opm-simulators/opm/simulators/wells/BlackoilWellModel_impl.hpp
Bård Skaflestad 6e9d2bd89e Initialise More Data Members at Construction Time
In particular, apply explicit default constructors to most data
members and push initialisation to initialiser list if convenient.

While here, also split long lines and apply const in more places.
Finally, reset well- and connection-level PI values to zero in
WellState::shutWell().  This is in preparation of including shut
wells in BlackoilWellModel's internal state.
2021-04-27 14:50:37 +02:00

3331 lines
128 KiB
C++

/*
Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics.
Copyright 2016 - 2018 Equinor ASA.
Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2016 - 2018 Norce AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <algorithm>
#include <utility>
#include <fmt/format.h>
namespace Opm {
template<typename TypeTag>
BlackoilWellModel<TypeTag>::
BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage)
: ebosSimulator_(ebosSimulator)
, terminal_output_((ebosSimulator.gridView().comm().rank() == 0) &&
EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput))
, phase_usage_(phase_usage)
, active_wgstate_(phase_usage)
, last_valid_wgstate_(phase_usage)
, nupcol_wgstate_(phase_usage)
{
// Create the guide rate container.
this->guideRate_ =
std::make_unique<GuideRate>(ebosSimulator_.vanguard().schedule());
local_num_cells_ = ebosSimulator_.gridView().size(0);
// Number of cells the global grid view
global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();
// Set up cartesian mapping.
{
const auto& grid = this->ebosSimulator_.vanguard().grid();
const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid),
cartDims[0] * cartDims[1] * cartDims[2]);
auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
this->parallel_well_info_.assign(parallel_wells.begin(),
parallel_wells.end());
}
const auto numProcs = ebosSimulator.gridView().comm().size();
this->not_on_process_ = [this, numProcs](const Well& well) {
if (well.getStatus() == Well::Status::SHUT)
return true;
if (numProcs == decltype(numProcs){1})
return false;
// Recall: false indicates NOT active!
const auto value = std::make_pair(well.name(), true);
auto candidate = std::lower_bound(this->parallel_well_info_.begin(),
this->parallel_well_info_.end(),
value);
return (candidate == this->parallel_well_info_.end())
|| (*candidate != value);
};
this->alternative_well_rate_init_ =
EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
}
template<typename TypeTag>
BlackoilWellModel<TypeTag>::
BlackoilWellModel(Simulator& ebosSimulator) :
BlackoilWellModel(ebosSimulator, phaseUsageFromDeck(ebosSimulator.vanguard().eclState()))
{}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
init()
{
extractLegacyCellPvtRegionIndex_();
extractLegacyDepth_();
gravity_ = ebosSimulator_.problem().gravity()[2];
initial_step_ = true;
// add the eWoms auxiliary module for the wells to the list
ebosSimulator_.model().addAuxiliaryModule(this);
is_cell_perforated_.resize(local_num_cells_, false);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
addNeighbors(std::vector<NeighborSet>& neighbors) const
{
if (!param_.matrix_add_well_contributions_) {
return;
}
// Create cartesian to compressed mapping
const auto& schedule_wells = schedule().getWellsatEnd();
// initialize the additional cell connections introduced by wells.
for (const auto& well : schedule_wells)
{
std::vector<int> wellCells;
// All possible connections of the well
const auto& connectionSet = well.getConnections();
wellCells.reserve(connectionSet.size());
for ( size_t c=0; c < connectionSet.size(); c++ )
{
const auto& connection = connectionSet.get(c);
int compressed_idx = cartesian_to_compressed_
.at(connection.global_index());
if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
wellCells.push_back(compressed_idx);
}
}
for (int cellIdx : wellCells) {
neighbors[cellIdx].insert(wellCells.begin(),
wellCells.end());
}
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
{
if (!localWellsActive())
return;
if (!param_.matrix_add_well_contributions_) {
// if the well contributions are not supposed to be included explicitly in
// the matrix, we only apply the vector part of the Schur complement here.
for (const auto& well: well_container_) {
// r = r - duneC_^T * invDuneD_ * resWell_
well->apply(res);
}
return;
}
for (const auto& well: well_container_) {
well->addWellContributions(jacobian);
// applying the well residual to reservoir residuals
// r = r - duneC_^T * invDuneD_ * resWell_
well->apply(res);
}
}
/// Return true if any well has a THP constraint.
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
hasTHPConstraints() const
{
int local_result = false;
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
for (const auto& well : well_container_) {
if (well->wellHasTHPConstraints(summaryState)) {
local_result=true;
}
}
return grid().comm().max(local_result);
}
/// Return true if the well was found and shut.
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
forceShutWellByNameIfPredictionMode(const std::string& wellname,
const double simulation_time)
{
// Only add the well to the closed list on the
// process that owns it.
int well_was_shut = 0;
for (const auto& well : well_container_) {
if (well->name() == wellname && !well->wellIsStopped()) {
if (well->underPredictionMode()) {
wellTestState_.closeWell(wellname, WellTestConfig::Reason::PHYSICAL, simulation_time);
well_was_shut = 1;
}
break;
}
}
// Communicate across processes if a well was shut.
well_was_shut = ebosSimulator_.vanguard().grid().comm().max(well_was_shut);
// Only log a message on the output rank.
if (terminal_output_ && well_was_shut) {
const std::string msg = "Well " + wellname
+ " will be shut because it cannot get converged.";
OpmLog::info(msg);
}
return (well_was_shut == 1);
}
template<typename TypeTag>
std::vector< Well >
BlackoilWellModel<TypeTag>::
getLocalWells(const int timeStepIdx, int& globalNumWells) const
{
auto w = schedule().getWells(timeStepIdx);
globalNumWells = w.size();
w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end());
return w;
}
template<typename TypeTag>
std::vector< ParallelWellInfo* >
BlackoilWellModel<TypeTag>::createLocalParallelWellInfo(const std::vector<Well>& wells)
{
std::vector< ParallelWellInfo* > local_parallel_well_info;
local_parallel_well_info.reserve(wells.size());
for (const auto& well : wells)
{
auto wellPair = std::make_pair(well.name(), true);
auto pwell = std::lower_bound(parallel_well_info_.begin(),
parallel_well_info_.end(),
wellPair);
assert(pwell != parallel_well_info_.end() &&
*pwell == wellPair);
local_parallel_well_info.push_back(&(*pwell));
}
return local_parallel_well_info;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
beginReportStep(const int timeStepIdx)
{
Opm::DeferredLogger local_deferredLogger;
report_step_starts_ = true;
const Grid& grid = ebosSimulator_.vanguard().grid();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
int globalNumWells = 0;
// Make wells_ecl_ contain only this partition's wells.
wells_ecl_ = getLocalWells(timeStepIdx, globalNumWells);
local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
// The well state initialize bhp with the cell pressure in the top cell.
// We must therefore provide it with updated cell pressures
this->initializeWellPerfData();
this->initializeWellState(timeStepIdx, globalNumWells, summaryState);
// Wells are active if they are active wells on at least
// one process.
wells_active_ = localWellsActive() ? 1 : 0;
wells_active_ = grid.comm().max(wells_active_);
// handling MS well related
if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
this->wellState().initWellStateMSWell(wells_ecl_, &this->prevWellState());
}
const Group& fieldGroup = schedule().getGroup("FIELD", timeStepIdx);
WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState());
// Compute reservoir volumes for RESV controls.
rateConverter_.reset(new RateConverterType (phase_usage_,
std::vector<int>(local_num_cells_, 0)));
rateConverter_->template defineState<ElementContext>(ebosSimulator_);
{
const auto& sched_state = this->schedule()[timeStepIdx];
// update VFP properties
vfp_properties_.reset(new VFPProperties( sched_state.vfpinj(), sched_state.vfpprod()) );
this->initializeWellProdIndCalculators();
if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) {
this->runWellPIScaling(timeStepIdx, local_deferredLogger);
}
}
// Store the current well state, to be able to recover in the case of failed iterations
this->commitWGState();
}
// called at the beginning of a time step
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
beginTimeStep()
{
updatePerforationIntensiveQuantities();
updateAverageFormationFactor();
Opm::DeferredLogger local_deferredLogger;
this->resetWGState();
updateAndCommunicateGroupData();
this->wellState().gliftTimeStepInit();
const int reportStepIdx = ebosSimulator_.episodeIndex();
const double simulationTime = ebosSimulator_.time();
std::string exc_msg;
auto exc_type = ExceptionType::NONE;
try {
// test wells
wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
// create the well container
well_container_ = createWellContainer(reportStepIdx);
// do the initialization for all the wells
// TODO: to see whether we can postpone of the intialization of the well containers to
// optimize the usage of the following several member variables
for (auto& well : well_container_) {
well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_);
}
// update the updated cell flag
std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
for (auto& well : well_container_) {
well->updatePerforatedCell(is_cell_perforated_);
}
// calculate the efficiency factors for each well
calculateEfficiencyFactors(reportStepIdx);
if constexpr (has_polymer_)
{
if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
setRepRadiusPerfLength();
}
}
} catch (const std::runtime_error& e) {
exc_type = ExceptionType::RUNTIME_ERROR;
exc_msg = e.what();
} catch (const std::invalid_argument& e) {
exc_type = ExceptionType::INVALID_ARGUMENT;
exc_msg = e.what();
} catch (const std::logic_error& e) {
exc_type = ExceptionType::LOGIC_ERROR;
exc_msg = e.what();
} catch (const std::exception& e) {
exc_type = ExceptionType::DEFAULT;
exc_msg = e.what();
}
logAndCheckForExceptionsAndThrow(local_deferredLogger, exc_type, "beginTimeStep() failed: " + exc_msg, terminal_output_);
for (auto& well : well_container_) {
well->setVFPProperties(vfp_properties_.get());
well->setGuideRate(guideRate_.get());
}
// Close completions due to economical reasons
for (auto& well : well_container_) {
well->closeCompletions(wellTestState_);
}
// calculate the well potentials
try {
std::vector<double> well_potentials;
computeWellPotentials(well_potentials, reportStepIdx, local_deferredLogger);
} catch ( std::runtime_error& e ) {
const std::string msg = "A zero well potential is returned for output purposes. ";
local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
}
if (alternative_well_rate_init_) {
// Update the well rates of well_state_, if only single-phase rates, to
// have proper multi-phase rates proportional to rates at bhp zero.
// This is done only for producers, as injectors will only have a single
// nonzero phase anyway.
for (auto& well : well_container_) {
if (well->isProducer()) {
well->updateWellStateRates(ebosSimulator_, this->wellState(), local_deferredLogger);
}
}
}
//update guide rates
const auto& comm = ebosSimulator_.vanguard().grid().comm();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
std::vector<double> pot(numPhases(), 0.0);
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
WellGroupHelpers::updateGuideRateForProductionGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, this->wellState(), this->groupState(), comm, guideRate_.get(), pot);
WellGroupHelpers::updateGuideRatesForInjectionGroups(fieldGroup, schedule(), summaryState, phase_usage_, reportStepIdx, this->wellState(), this->groupState(), guideRate_.get(), local_deferredLogger);
WellGroupHelpers::updateGuideRatesForWells(schedule(), phase_usage_, reportStepIdx, simulationTime, this->wellState(), comm, guideRate_.get());
try {
// Compute initial well solution for new wells and injectors that change injection type i.e. WAG.
for (auto& well : well_container_) {
const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+ ScheduleEvents::INJECTION_TYPE_CHANGED
+ ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
+ ScheduleEvents::NEW_WELL;
const auto& events = schedule()[reportStepIdx].wellgroup_events();
const bool event = report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
if (event) {
try {
well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), local_deferredLogger);
well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), local_deferredLogger);
} catch (const std::exception& e) {
const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates";
local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
}
}
}
} catch (const std::runtime_error& e) {
exc_type = ExceptionType::RUNTIME_ERROR;
exc_msg = e.what();
} catch (const std::invalid_argument& e) {
exc_type = ExceptionType::INVALID_ARGUMENT;
exc_msg = e.what();
} catch (const std::logic_error& e) {
exc_type = ExceptionType::LOGIC_ERROR;
exc_msg = e.what();
} catch (const std::exception& e) {
exc_type = ExceptionType::DEFAULT;
exc_msg = e.what();
}
if (exc_type != ExceptionType::NONE) {
const std::string msg = "Compute initial well solution for new wells failed. Continue with zero initial rates";
local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
}
logAndCheckForExceptionsAndThrow(local_deferredLogger,
exc_type, "beginTimeStep() failed: " + exc_msg, terminal_output_);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::gliftDebug(
const std::string &msg, Opm::DeferredLogger &deferred_logger) const
{
if (this->glift_debug) {
const std::string message = fmt::format(
" GLIFT (DEBUG) : BlackoilWellModel : {}", msg);
deferred_logger.info(message);
}
}
template<typename TypeTag>
void
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
const auto wellsForTesting = wellTestState_
.updateWells(wtest_config, wells_ecl_, simulationTime);
for (const auto& testWell : wellsForTesting) {
const std::string& well_name = testWell.first;
// this is the well we will test
WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
// 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);
well->setWellEfficiencyFactor(well_efficiency_factor);
well->setVFPProperties(vfp_properties_.get());
well->setGuideRate(guideRate_.get());
const WellTestConfig::Reason testing_reason = testWell.second;
well->wellTesting(ebosSimulator_, simulationTime, timeStepIdx,
testing_reason, this->wellState(), this->groupState(), wellTestState_, deferred_logger);
}
}
}
// called at the end of a report step
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
endReportStep()
{
// Clear the communication data structures for above values.
for (auto&& pinfo : local_parallel_well_info_)
{
pinfo->clear();
}
}
// 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)
{
this->closed_this_step_.clear();
// time step is finished and we are not any more at the beginning of an report step
report_step_starts_ = false;
const int reportStepIdx = ebosSimulator_.episodeIndex();
Opm::DeferredLogger local_deferredLogger;
for (const auto& well : well_container_) {
if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
well->updateWaterThroughput(dt, this->wellState());
}
}
updateWellTestState(simulationTime, wellTestState_);
// update the rate converter with current averages pressures etc in
rateConverter_->template defineState<ElementContext>(ebosSimulator_);
// calculate the well potentials
try {
std::vector<double> well_potentials;
computeWellPotentials(well_potentials, reportStepIdx, local_deferredLogger);
} catch ( std::runtime_error& e ) {
const std::string msg = "A zero well potential is returned for output purposes. ";
local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
}
// check group sales limits at the end of the timestep
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
checkGconsaleLimits(fieldGroup, this->wellState(), local_deferredLogger);
this->calculateProductivityIndexValues(local_deferredLogger);
this->commitWGState();
Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
if (terminal_output_) {
global_deferredLogger.logMessages();
}
//reporting output temperatures
this->computeWellTemperature();
}
template<typename TypeTag>
template <class Context>
void
BlackoilWellModel<TypeTag>::
computeTotalRatesForDof(RateVector& rate,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx) const
{
rate = 0;
int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
if (!is_cell_perforated_[elemIdx])
return;
for (const auto& well : well_container_)
well->addCellRates(rate, elemIdx);
}
template<typename TypeTag>
typename BlackoilWellModel<TypeTag>::WellInterfacePtr
BlackoilWellModel<TypeTag>::
well(const std::string& wellName) const
{
for (const auto& well : well_container_) {
if (well->name() == wellName) {
return well;
}
}
OPM_THROW(std::invalid_argument, "The well with name " + wellName + " is not in the well Container");
return nullptr;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initFromRestartFile(const RestartValue& restartValues)
{
// The restart step value is used to identify wells present at the given
// time step. Wells that are added at the same time step as RESTART is initiated
// will not be present in a restart file. Use the previous time step to retrieve
// wells that have information written to the restart file.
const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0);
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
int globalNumWells = 0;
// wells_ecl_ should only contain wells on this processor.
wells_ecl_ = getLocalWells(report_step, globalNumWells);
local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
this->initializeWellProdIndCalculators();
initializeWellPerfData();
const int nw = wells_ecl_.size();
if (nw > 0) {
const auto phaseUsage = phaseUsageFromDeck(eclState());
const size_t numCells = Opm::UgGridHelpers::numCells(grid());
const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal());
this->wellState().resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, well_perf_data_, summaryState, globalNumWells); // Resize for restart step
loadRestartData(restartValues.wells, restartValues.grp_nwrk, phaseUsage, handle_ms_well, this->wellState());
}
this->commitWGState();
initial_step_ = false;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initializeWellProdIndCalculators()
{
this->prod_index_calc_.clear();
this->prod_index_calc_.reserve(this->wells_ecl_.size());
for (const auto& well : this->wells_ecl_) {
this->prod_index_calc_.emplace_back(well);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initializeWellPerfData()
{
well_perf_data_.resize(wells_ecl_.size());
int well_index = 0;
for (const auto& well : wells_ecl_) {
int completion_index = 0;
// INVALID_ECL_INDEX marks no above perf available
int completion_index_above = ParallelWellInfo::INVALID_ECL_INDEX;
well_perf_data_[well_index].clear();
well_perf_data_[well_index].reserve(well.getConnections().size());
CheckDistributedWellConnections checker(well, *local_parallel_well_info_[well_index]);
bool hasFirstPerforation = false;
bool firstOpenCompletion = true;
auto& parallelWellInfo = *local_parallel_well_info_[well_index];
parallelWellInfo.beginReset();
for (const auto& completion : well.getConnections()) {
const int active_index =
cartesian_to_compressed_[completion.global_index()];
if (completion.state() == Connection::State::OPEN) {
if (active_index >= 0) {
if (firstOpenCompletion)
{
hasFirstPerforation = true;
}
checker.connectionFound(completion_index);
PerforationData pd;
pd.cell_index = active_index;
pd.connection_transmissibility_factor = completion.CF();
pd.satnum_id = completion.satTableId();
pd.ecl_index = completion_index;
well_perf_data_[well_index].push_back(pd);
parallelWellInfo.pushBackEclIndex(completion_index_above,
completion_index);
}
firstOpenCompletion = false;
// Next time this index is the one above as each open completion is
// is stored somehwere.
completion_index_above = completion_index;
} else {
checker.connectionFound(completion_index);
if (completion.state() != Connection::State::SHUT) {
OPM_THROW(std::runtime_error,
"Completion state: " << Connection::State2String(completion.state()) << " not handled");
}
}
// Note: we rely on the connections being filtered! I.e. there are only connections
// to active cells in the global grid.
++completion_index;
}
parallelWellInfo.endReset();
checker.checkAllConnectionsFound();
parallelWellInfo.communicateFirstPerforation(hasFirstPerforation);
++well_index;
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initializeWellState(const int timeStepIdx,
const int globalNumWells,
const SummaryState& summaryState)
{
std::vector<double> cellPressures(this->local_num_cells_, 0.0);
ElementContext elemCtx(ebosSimulator_);
const auto& gridView = ebosSimulator_.vanguard().gridView();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
if (elemIt->partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState();
// copy of get perfpressure in Standard well except for value
double& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)];
if (Indices::oilEnabled) {
perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
} else if (Indices::waterEnabled) {
perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
} else {
perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
}
}
this->wellState().init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx,
&this->prevWellState(), well_perf_data_,
summaryState, globalNumWells);
}
template<typename TypeTag>
std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
BlackoilWellModel<TypeTag>::
createWellContainer(const int time_step)
{
std::vector<WellInterfacePtr> well_container;
Opm::DeferredLogger local_deferredLogger;
const int nw = numLocalWells();
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();
const auto well_status = this->schedule()
.getWell(well_name, time_step).getStatus();
if ((well_ecl.getStatus() == Well::Status::SHUT) ||
(well_status == Well::Status::SHUT))
{
// Due to ACTIONX the well might have been closed behind our back.
if (well_ecl.getStatus() != Well::Status::SHUT) {
this->closed_this_step_.insert(well_name);
this->wellState().shutWell(w);
}
continue;
}
// A new WCON keywords can re-open a well that was closed/shut due to Physical limit
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 (this->wellState().effectiveEventsOccurred(w)) {
if (wellTestState_.lastTestTime(well_name) == ebosSimulator_.time()) {
// The well was shut this timestep, we are most likely retrying
// a timestep without the well in question, after it caused
// repeated timestep cuts. It should therefore not be opened,
// even if it was new or received new targets this report step.
this->wellState().setEffectiveEventsOccurred(w, false);
} else {
wellTestState_.openWell(well_name);
}
}
}
// TODO: should we do this for all kinds of closing reasons?
// something like wellTestState_.hasWell(well_name)?
bool wellIsStopped = false;
if (wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::ECONOMIC) ||
wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::PHYSICAL))
{
if (well_ecl.getAutomaticShutIn()) {
// shut wells are not added to the well container
this->wellState().shutWell(w);
continue;
} else {
// stopped wells are added to the container but marked as stopped
this->wellState().stopWell(w);
wellIsStopped = true;
}
}
// If a production well disallows crossflow and its
// (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();
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 = is_zero(prod_controls.oil_rate);
break;
case Well::ProducerCMode::WRAT:
zero_rate_control = is_zero(prod_controls.water_rate);
break;
case Well::ProducerCMode::GRAT:
zero_rate_control = is_zero(prod_controls.gas_rate);
break;
case Well::ProducerCMode::LRAT:
zero_rate_control = is_zero(prod_controls.liquid_rate);
break;
case Well::ProducerCMode::RESV:
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());
this->wellState().shutWell(w);
continue;
}
}
if (well_status == Well::Status::STOP) {
this->wellState().stopWell(w);
wellIsStopped = true;
}
well_container.emplace_back(this->createWellPointer(w, time_step));
if (wellIsStopped)
well_container.back()->stopWell();
}
}
// Collect log messages and print.
Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
if (terminal_output_) {
global_deferredLogger.logMessages();
}
return well_container;
}
template <typename TypeTag>
typename BlackoilWellModel<TypeTag>::WellInterfacePtr
BlackoilWellModel<TypeTag>::
createWellPointer(const int wellID, const int time_step) const
{
const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
return this->template createTypedWellPointer<StandardWell<TypeTag>>(wellID, time_step);
}
else {
return this->template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, time_step);
}
}
template <typename TypeTag>
template <typename WellType>
std::unique_ptr<WellType>
BlackoilWellModel<TypeTag>::
createTypedWellPointer(const int wellID, const int time_step) const
{
// Use the pvtRegionIdx from the top cell
const auto& perf_data = this->well_perf_data_[wellID];
// Cater for case where local part might have no perforations.
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];
const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
return std::make_unique<WellType>(this->wells_ecl_[wellID],
parallel_well_info,
time_step,
this->param_,
*this->rateConverter_,
global_pvtreg,
this->numComponents(),
this->numPhases(),
wellID,
this->wellState().firstPerfIndex()[wellID],
perf_data);
}
template<typename TypeTag>
typename BlackoilWellModel<TypeTag>::WellInterfacePtr
BlackoilWellModel<TypeTag>::
createWellForWellTest(const std::string& well_name,
const int report_step,
Opm::DeferredLogger& deferred_logger) const
{
// Finding the location of the well in wells_ecl
const int nw_wells_ecl = wells_ecl_.size();
int index_well_ecl = 0;
for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
if (well_name == wells_ecl_[index_well_ecl].name()) {
break;
}
}
// It should be able to find in wells_ecl.
if (index_well_ecl == nw_wells_ecl) {
OPM_DEFLOG_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ", deferred_logger);
}
return this->createWellPointer(index_well_ecl, report_step);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assemble(const int iterationIdx,
const double dt)
{
Opm::DeferredLogger local_deferredLogger;
if (this->glift_debug) {
const std::string msg = fmt::format(
"assemble() : iteration {}" , iterationIdx);
gliftDebug(msg, local_deferredLogger);
}
last_report_ = SimulatorReportSingle();
Dune::Timer perfTimer;
perfTimer.start();
if ( ! wellsActive() ) {
return;
}
updatePerforationIntensiveQuantities();
auto exc_type = ExceptionType::NONE;
std::string exc_msg;
try {
if (iterationIdx == 0) {
calculateExplicitQuantities(local_deferredLogger);
prepareTimeStep(local_deferredLogger);
}
updateWellControls(local_deferredLogger, /* check group controls */ true);
// Set the well primary variables based on the value of well solutions
initPrimaryVariablesEvaluation();
if (param_.solve_welleq_initially_ && iterationIdx == 0) {
for (auto& well : well_container_) {
well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), local_deferredLogger);
}
updateWellControls(local_deferredLogger, /* check group controls */ false);
}
maybeDoGasLiftOptimize(local_deferredLogger);
assembleWellEq(dt, local_deferredLogger);
} catch (const std::runtime_error& e) {
exc_type = ExceptionType::RUNTIME_ERROR;
exc_msg = e.what();
} catch (const std::invalid_argument& e) {
exc_type = ExceptionType::INVALID_ARGUMENT;
exc_msg = e.what();
} catch (const std::logic_error& e) {
exc_type = ExceptionType::LOGIC_ERROR;
exc_msg = e.what();
} catch (const std::exception& e) {
exc_type = ExceptionType::DEFAULT;
exc_msg = e.what();
}
logAndCheckForExceptionsAndThrow(local_deferredLogger, exc_type, "assemble() failed: " + exc_msg, terminal_output_);
last_report_.converged = true;
last_report_.assemble_time_well += perfTimer.stop();
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
maybeDoGasLiftOptimize(Opm::DeferredLogger& deferred_logger)
{
this->wellState().enableGliftOptimization();
GLiftOptWells glift_wells;
GLiftProdWells prod_wells;
GLiftWellStateMap state_map;
// Stage1: Optimize single wells not checking any group limits
for (auto& well : well_container_) {
well->gasLiftOptimizationStage1(
this->wellState(), ebosSimulator_, deferred_logger,
prod_wells, glift_wells, state_map);
}
gasLiftOptimizationStage2(deferred_logger, prod_wells, glift_wells, state_map);
if (this->glift_debug) gliftDebugShowALQ(deferred_logger);
this->wellState().disableGliftOptimization();
}
// If a group has any production rate constraints, and/or a limit
// on its total rate of lift gas supply, allocate lift gas
// preferentially to the wells that gain the most benefit from
// it. Lift gas increments are allocated in turn to the well that
// currently has the largest weighted incremental gradient. The
// procedure takes account of any limits on the group production
// rate or lift gas supply applied to any level of group.
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
gasLiftOptimizationStage2(Opm::DeferredLogger& deferred_logger,
GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
GLiftWellStateMap &glift_well_state_map)
{
GasLiftStage2 glift {*this, ebosSimulator_, deferred_logger, this->wellState(),
prod_wells, glift_wells, glift_well_state_map};
glift.runOptimize();
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
gliftDebugShowALQ(Opm::DeferredLogger& deferred_logger)
{
for (auto& well : this->well_container_) {
if (well->isProducer()) {
auto alq = this->wellState().getALQ(well->name());
const std::string msg = fmt::format("ALQ_REPORT : {} : {}",
well->name(), alq);
gliftDebug(msg, deferred_logger);
}
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assembleWellEq(const double dt, Opm::DeferredLogger& deferred_logger)
{
for (auto& well : well_container_) {
well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
apply( BVector& r) const
{
if ( ! localWellsActive() ) {
return;
}
for (auto& well : well_container_) {
well->apply(r);
}
}
// Ax = A x - C D^-1 B x
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
apply(const BVector& x, BVector& Ax) const
{
// TODO: do we still need localWellsActive()?
if ( ! localWellsActive() ) {
return;
}
for (auto& well : well_container_) {
well->apply(x, Ax);
}
}
#if HAVE_CUDA || HAVE_OPENCL
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
getWellContributions(WellContributions& wellContribs) const
{
// prepare for StandardWells
wellContribs.setBlockSize(StandardWell<TypeTag>::numEq, StandardWell<TypeTag>::numStaticWellEq);
for(unsigned int i = 0; i < well_container_.size(); i++){
auto& well = well_container_[i];
std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
if (derived) {
unsigned int numBlocks;
derived->getNumBlocks(numBlocks);
wellContribs.addNumBlocks(numBlocks);
}
}
// allocate memory for data from StandardWells
wellContribs.alloc();
for(unsigned int i = 0; i < well_container_.size(); i++){
auto& well = well_container_[i];
// maybe WellInterface could implement addWellContribution()
auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
if (derived_std) {
derived_std->addWellContribution(wellContribs);
} else {
auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
if (derived_ms) {
derived_ms->addWellContribution(wellContribs);
} else {
OpmLog::warning("Warning unknown type of well");
}
}
}
}
#endif
// Ax = Ax - alpha * C D^-1 B x
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
{
if ( ! localWellsActive() ) {
return;
}
if( scaleAddRes_.size() != Ax.size() ) {
scaleAddRes_.resize( Ax.size() );
}
scaleAddRes_ = 0.0;
// scaleAddRes_ = - C D^-1 B x
apply( x, scaleAddRes_ );
// Ax = Ax + alpha * scaleAddRes_
Ax.axpy( alpha, scaleAddRes_ );
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
recoverWellSolutionAndUpdateWellState(const BVector& x)
{
Opm::DeferredLogger local_deferredLogger;
auto exc_type = ExceptionType::NONE;
std::string exc_msg;
try {
if (localWellsActive()) {
for (auto& well : well_container_) {
well->recoverWellSolutionAndUpdateWellState(x, this->wellState(), local_deferredLogger);
}
}
} catch (const std::runtime_error& e) {
exc_type = ExceptionType::RUNTIME_ERROR;
exc_msg = e.what();
} catch (const std::invalid_argument& e) {
exc_type = ExceptionType::INVALID_ARGUMENT;
exc_msg = e.what();
} catch (const std::logic_error& e) {
exc_type = ExceptionType::LOGIC_ERROR;
exc_msg = e.what();
} catch (const std::exception& e) {
exc_type = ExceptionType::DEFAULT;
exc_msg = e.what();
}
logAndCheckForExceptionsAndThrow(local_deferredLogger, exc_type, "recoverWellSolutionAndUpdateWellState() failed: " + exc_msg, terminal_output_);
}
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
wellsActive() const
{
return wells_active_;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
setWellsActive(const bool wells_active)
{
wells_active_ = wells_active;
}
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
localWellsActive() const
{
return numLocalWells() > 0;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initPrimaryVariablesEvaluation() const
{
for (auto& well : well_container_) {
well->initPrimaryVariablesEvaluation();
}
}
template<typename TypeTag>
ConvergenceReport
BlackoilWellModel<TypeTag>::
getWellConvergence(const std::vector<Scalar>& B_avg, bool checkGroupConvergence) const
{
Opm::DeferredLogger local_deferredLogger;
// Get global (from all processes) convergence report.
ConvergenceReport local_report;
for (const auto& well : well_container_) {
if (well->isOperable() ) {
local_report += well->getWellConvergence(this->wellState(), B_avg, local_deferredLogger);
}
}
Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
if (terminal_output_) {
global_deferredLogger.logMessages();
}
ConvergenceReport report = gatherConvergenceReport(local_report);
// Log debug messages for NaN or too large residuals.
if (terminal_output_) {
for (const auto& f : report.wellFailures()) {
if (f.severity() == ConvergenceReport::Severity::NotANumber) {
OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
} else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
}
}
}
if (checkGroupConvergence) {
const int reportStepIdx = ebosSimulator_.episodeIndex();
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
bool violated = checkGroupConstraints(fieldGroup, global_deferredLogger);
report.setGroupConverged(!violated);
}
return report;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
calculateExplicitQuantities(Opm::DeferredLogger& deferred_logger) const
{
// TODO: checking isOperable() ?
for (auto& well : well_container_) {
well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateWellControls(Opm::DeferredLogger& deferred_logger, const bool checkGroupControls)
{
// Even if there are no wells active locally, we cannot
// return as the DeferredLogger uses global communication.
// For no well active globally we simply return.
if( !wellsActive() ) return ;
updateAndCommunicateGroupData();
updateNetworkPressures();
std::set<std::string> switched_wells;
std::set<std::string> switched_groups;
if (checkGroupControls) {
// Check group individual constraints.
updateGroupIndividualControls(deferred_logger, switched_groups);
// Check group's constraints from higher levels.
updateGroupHigherControls(deferred_logger, switched_groups);
updateAndCommunicateGroupData();
// Check wells' group constraints and communicate.
for (const auto& well : well_container_) {
const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
const bool changed = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
if (changed) {
switched_wells.insert(well->name());
}
}
updateAndCommunicateGroupData();
}
// Check individual well constraints and communicate.
for (const auto& well : well_container_) {
if (switched_wells.count(well->name())) {
continue;
}
const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
}
updateAndCommunicateGroupData();
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateNetworkPressures()
{
// Get the network and return if inactive.
const int reportStepIdx = ebosSimulator_.episodeIndex();
const auto& network = schedule()[reportStepIdx].network();
if (!network.active()) {
return;
}
node_pressures_ = WellGroupHelpers::computeNetworkPressures(network, this->wellState(), this->groupState(), *(vfp_properties_->getProd()), schedule(), reportStepIdx);
// Set the thp limits of wells
for (auto& well : well_container_) {
// Producers only, since we so far only support the
// "extended" network model (properties defined by
// BRANPROP and NODEPROP) which only applies to producers.
if (well->isProducer()) {
const auto it = node_pressures_.find(well->wellEcl().groupName());
if (it != node_pressures_.end()) {
// The well belongs to a group with has a network pressure constraint,
// set the dynamic THP constraint of the well accordingly.
well->setDynamicThpLimit(it->second);
}
}
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateAndCommunicateGroupData()
{
const int reportStepIdx = ebosSimulator_.episodeIndex();
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
const int nupcol = schedule()[reportStepIdx].nupcol();
const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
// This builds some necessary lookup structures, so it must be called
// before we copy to well_state_nupcol_.
const auto& comm = ebosSimulator_.vanguard().grid().comm();
this->wellState().updateGlobalIsGrup(schedule(), reportStepIdx, comm);
if (iterationIdx < nupcol) {
this->updateNupcolWGState();
}
auto& well_state = this->wellState();
const auto& well_state_nupcol = this->nupcolWellState();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
// the group target reduction rates needs to be update since wells may have swicthed to/from GRUP control
// Currently the group target reduction does not honor NUPCOL. TODO: is that true?
std::vector<double> groupTargetReduction(numPhases(), 0.0);
WellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, phase_usage_, *guideRate_, well_state_nupcol, well_state, this->groupState(), groupTargetReduction);
std::vector<double> groupTargetReductionInj(numPhases(), 0.0);
WellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, phase_usage_, *guideRate_, well_state_nupcol, well_state, this->groupState(), groupTargetReductionInj);
WellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, phase_usage_, summaryState, well_state_nupcol, well_state, this->groupState());
WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState());
WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState());
WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState());
// We use the rates from the privious time-step to reduce oscilations
WellGroupHelpers::updateWellRates(fieldGroup, schedule(), reportStepIdx, this->prevWellState(), well_state);
// Set ALQ for off-process wells to zero
for (const auto& wname : schedule().wellNames(reportStepIdx)) {
const bool is_producer = schedule().getWell(wname, reportStepIdx).isProducer();
const bool not_on_this_process = well_state.wellMap().count(wname) == 0;
if (is_producer && not_on_this_process) {
well_state.setALQ(wname, 0.0);
}
}
well_state.communicateGroupRates(comm);
this->groupState().communicate_rates(comm);
// compute wsolvent fraction for REIN wells
updateWsolvent(fieldGroup, schedule(), reportStepIdx, well_state_nupcol);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
{
Opm::DeferredLogger local_deferredLogger;
for (const auto& well : well_container_) {
const auto wasClosed = wellTestState.hasWellClosed(well->name());
well->updateWellTestState(this->wellState(), simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
if (!wasClosed && wellTestState.hasWellClosed(well->name())) {
this->closed_this_step_.insert(well->name());
}
}
Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
if (terminal_output_) {
global_deferredLogger.logMessages();
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
computeWellPotentials(std::vector<double>& well_potentials, const int reportStepIdx, Opm::DeferredLogger& deferred_logger)
{
// number of wells and phases
const int nw = numLocalWells();
const int np = numPhases();
well_potentials.resize(nw * np, 0.0);
auto well_state = this->wellState();
const Opm::SummaryConfig& summaryConfig = ebosSimulator_.vanguard().summaryConfig();
const bool write_restart_file = ebosSimulator_.vanguard().schedule().write_rst_file(reportStepIdx);
auto exc_type = ExceptionType::NONE;
std::string exc_msg;
for (const auto& well : well_container_) {
const bool needed_for_summary = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) ||
summaryConfig.hasSummaryKey( "WOPI:" + well->name()) ||
summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->isInjector()) ||
((summaryConfig.hasSummaryKey( "WWPP:" + well->name()) ||
summaryConfig.hasSummaryKey( "WOPP:" + well->name()) ||
summaryConfig.hasSummaryKey( "WGPP:" + well->name())) && well->isProducer());
bool needPotentialsForGuideRate = true;//eclWell.getGuideRatePhase() == Well::GuideRateTarget::UNDEFINED;
if (write_restart_file || needed_for_summary || needPotentialsForGuideRate)
{
try {
std::vector<double> potentials;
well->computeWellPotentials(ebosSimulator_, well_state, potentials, deferred_logger);
// putting the sucessfully calculated potentials to the well_potentials
for (int p = 0; p < np; ++p) {
well_potentials[well->indexOfWell() * np + p] = std::abs(potentials[p]);
}
} catch (const std::runtime_error& e) {
exc_type = ExceptionType::RUNTIME_ERROR;
exc_msg = e.what();
} catch (const std::invalid_argument& e) {
exc_type = ExceptionType::INVALID_ARGUMENT;
exc_msg = e.what();
} catch (const std::logic_error& e) {
exc_type = ExceptionType::LOGIC_ERROR;
exc_msg = e.what();
} catch (const std::exception& e) {
exc_type = ExceptionType::DEFAULT;
exc_msg = e.what();
}
}
}
logAndCheckForExceptionsAndThrow(deferred_logger, exc_type,
"computeWellPotentials() failed: " + exc_msg,
terminal_output_);
// Store it in the well state
this->wellState().wellPotentials() = well_potentials;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
calculateProductivityIndexValues(DeferredLogger& deferred_logger)
{
if (! this->localWellsActive()) {
return;
}
for (const auto& wellPtr : this->well_container_) {
wellPtr->updateProductivityIndex(this->ebosSimulator_,
this->prod_index_calc_[wellPtr->indexOfWell()],
this->wellState(),
deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
prepareTimeStep(Opm::DeferredLogger& deferred_logger)
{
auto exc_type = ExceptionType::NONE;
std::string exc_msg;
try {
for (const auto& well : well_container_) {
well->checkWellOperability(ebosSimulator_, this->wellState(), deferred_logger);
}
// since the controls are all updated, we should update well_state accordingly
for (const auto& well : well_container_) {
const int w = well->indexOfWell();
if (!well->isOperable() ) continue;
if (this->wellState().effectiveEventsOccurred(w) ) {
well->updateWellStateWithTarget(ebosSimulator_, this->wellState(), deferred_logger);
}
// there is no new well control change input within a report step,
// so next time step, the well does not consider to have effective events anymore
// TODO: if we can know whether this is the first time step within the report step,
// we do not need to set it to false
// TODO: we should do this at the end of the time step in case we will need it within
// this time step somewhere
if (this->wellState().effectiveEventsOccurred(w) ) {
this->wellState().setEffectiveEventsOccurred(w, false);
}
} // end of for (const auto& well : well_container_)
updatePrimaryVariables(deferred_logger);
} catch (const std::runtime_error& e) {
exc_type = ExceptionType::RUNTIME_ERROR;
exc_msg = e.what();
} catch (const std::invalid_argument& e) {
exc_type = ExceptionType::INVALID_ARGUMENT;
exc_msg = e.what();
} catch (const std::logic_error& e) {
exc_type = ExceptionType::LOGIC_ERROR;
exc_msg = e.what();
} catch (const std::exception& e) {
exc_type = ExceptionType::DEFAULT;
exc_msg = e.what();
}
logAndCheckForExceptionsAndThrow(deferred_logger, exc_type, "prepareTimestep() failed: " + exc_msg, terminal_output_);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
calculateEfficiencyFactors(const int reportStepIdx)
{
if ( !localWellsActive() ) {
return;
}
for (auto& well : well_container_) {
const Well& wellEcl = well->wellEcl();
double well_efficiency_factor = wellEcl.getEfficiencyFactor();
WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), reportStepIdx), schedule(), reportStepIdx, well_efficiency_factor);
well->setWellEfficiencyFactor(well_efficiency_factor);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
setupCartesianToCompressed_(const int* global_cell, int number_of_cartesian_cells)
{
cartesian_to_compressed_.resize(number_of_cartesian_cells, -1);
if (global_cell) {
auto elemIt = ebosSimulator_.gridView().template begin</*codim=*/ 0>();
for (unsigned i = 0; i < local_num_cells_; ++i) {
// Skip perforations in the overlap/ghost for distributed wells.
if (elemIt->partitionType() == Dune::InteriorEntity)
{
assert(ebosSimulator_.gridView().indexSet().index(*elemIt) == static_cast<int>(i));
cartesian_to_compressed_[global_cell[i]] = i;
}
++elemIt;
}
}
else {
for (unsigned i = 0; i < local_num_cells_; ++i) {
cartesian_to_compressed_[i] = i;
}
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
setRepRadiusPerfLength()
{
for (const auto& well : well_container_) {
well->setRepRadiusPerfLength(cartesian_to_compressed_);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateAverageFormationFactor()
{
std::vector< Scalar > B_avg(numComponents(), Scalar() );
const auto& grid = ebosSimulator_.vanguard().grid();
const auto& gridView = grid.leafGridView();
ElementContext elemCtx(ebosSimulator_);
const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
elemIt != elemEndIt; ++elemIt)
{
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
{
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
auto& B = B_avg[ compIdx ];
B += 1 / fs.invB(phaseIdx).value();
}
if constexpr (has_solvent_) {
auto& B = B_avg[solventSaturationIdx];
B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
}
}
// compute global average
grid.comm().sum(B_avg.data(), B_avg.size());
for(auto& bval: B_avg)
{
bval/=global_num_cells_;
}
B_avg_ = B_avg;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updatePrimaryVariables(Opm::DeferredLogger& deferred_logger)
{
for (const auto& well : well_container_) {
well->updatePrimaryVariables(this->wellState(), deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::extractLegacyCellPvtRegionIndex_()
{
const auto& grid = ebosSimulator_.vanguard().grid();
const auto& eclProblem = ebosSimulator_.problem();
const unsigned numCells = grid.size(/*codim=*/0);
pvt_region_idx_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
pvt_region_idx_[cellIdx] =
eclProblem.pvtRegionIndex(cellIdx);
}
}
// The number of components in the model.
template<typename TypeTag>
int
BlackoilWellModel<TypeTag>::numComponents() const
{
if (wellsActive() && numPhases() < 3) {
return numPhases();
}
int numComp = FluidSystem::numComponents;
if constexpr (has_solvent_) {
numComp ++;
}
return numComp;
}
template<typename TypeTag>
int
BlackoilWellModel<TypeTag>:: numLocalWells() const
{
return wells_ecl_.size();
}
template<typename TypeTag>
int
BlackoilWellModel<TypeTag>::numPhases() const
{
return phase_usage_.num_phases;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::extractLegacyDepth_()
{
const auto& grid = ebosSimulator_.vanguard().grid();
const unsigned numCells = grid.size(/*codim=*/0);
depth_.resize(numCells);
for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
depth_[cellIdx] = Opm::UgGridHelpers::cellCenterDepth( grid, cellIdx );
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updatePerforationIntensiveQuantities() {
ElementContext elemCtx(ebosSimulator_);
const auto& gridView = ebosSimulator_.gridView();
const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
elemIt != elemEndIt;
++elemIt)
{
elemCtx.updatePrimaryStencil(*elemIt);
int elemIdx = elemCtx.globalSpaceIndex(0, 0);
if (!is_cell_perforated_[elemIdx]) {
continue;
}
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
}
}
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
hasWell(const std::string& wname) {
auto iter = std::find_if( this->wells_ecl_.begin(), this->wells_ecl_.end(), [&wname](const Well& well) { return well.name() == wname; });
return (iter != this->wells_ecl_.end());
}
// convert well data from opm-common to well state from opm-core
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
loadRestartData( const data::Wells& rst_wells,
const data::GroupAndNetworkValues& grpNwrkValues,
const PhaseUsage& phases,
const bool handle_ms_well,
WellStateFullyImplicitBlackoil& well_state)
{
using GPMode = Group::ProductionCMode;
using GIMode = Group::InjectionCMode;
using rt = data::Rates::opt;
const auto np = phases.num_phases;
std::vector< rt > phs( np );
if( phases.phase_used[BlackoilPhases::Aqua] ) {
phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat;
}
if( phases.phase_used[BlackoilPhases::Liquid] ) {
phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil;
}
if( phases.phase_used[BlackoilPhases::Vapour] ) {
phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
}
for( const auto& wm : well_state.wellMap() ) {
const auto well_index = wm.second[ 0 ];
const auto& rst_well = rst_wells.at( wm.first );
well_state.bhp()[ well_index ] = rst_well.bhp;
well_state.temperature()[ well_index ] = rst_well.temperature;
if (rst_well.current_control.isProducer) {
well_state.currentProductionControls()[ well_index ] = rst_well.current_control.prod;
}
else {
well_state.currentInjectionControls()[ well_index ] = rst_well.current_control.inj;
}
const auto wellrate_index = well_index * np;
for( size_t i = 0; i < phs.size(); ++i ) {
assert( rst_well.rates.has( phs[ i ] ) );
well_state.wellRates()[ wellrate_index + i ] = rst_well.rates.get( phs[ i ] );
}
auto * perf_pressure = well_state.perfPress().data() + wm.second[1];
auto * perf_rates = well_state.perfRates().data() + wm.second[1];
auto * perf_phase_rates = well_state.perfPhaseRates().data() + wm.second[1]*np;
const auto& perf_data = this->well_perf_data_[well_index];
for (std::size_t perf_index = 0; perf_index < perf_data.size(); perf_index++) {
const auto& pd = perf_data[perf_index];
const auto& rst_connection = rst_well.connections[pd.ecl_index];
perf_pressure[perf_index] = rst_connection.pressure;
perf_rates[perf_index] = rst_connection.reservoir_rate;
for (int phase_index = 0; phase_index < np; ++phase_index)
perf_phase_rates[perf_index*np + phase_index] = rst_connection.rates.get(phs[phase_index]);
}
if (handle_ms_well && !rst_well.segments.empty()) {
// we need the well_ecl_ information
const std::string& well_name = wm.first;
const Well& well_ecl = getWellEcl(well_name);
const WellSegments& segment_set = well_ecl.getSegments();
const int top_segment_index = well_state.topSegmentIndex(well_index);
const auto& segments = rst_well.segments;
// \Note: eventually we need to hanlde the situations that some segments are shut
assert(0u + segment_set.size() == segments.size());
for (const auto& segment : segments) {
const int segment_index = segment_set.segmentNumberToIndex(segment.first);
// recovering segment rates and pressure from the restart values
const auto pres_idx = Opm::data::SegmentPressures::Value::Pressure;
well_state.segPress()[top_segment_index + segment_index] = segment.second.pressures[pres_idx];
const auto& segment_rates = segment.second.rates;
for (int p = 0; p < np; ++p) {
well_state.segRates()[(top_segment_index + segment_index) * np + p] = segment_rates.get(phs[p]);
}
}
}
}
for (const auto& [group, value] : grpNwrkValues.groupData) {
const auto cpc = value.currentControl.currentProdConstraint;
const auto cgi = value.currentControl.currentGasInjectionConstraint;
const auto cwi = value.currentControl.currentWaterInjectionConstraint;
if (cpc != GPMode::NONE) {
this->groupState().production_control(group, cpc);
}
if (cgi != GIMode::NONE) {
this->groupState().injection_control(group, Phase::GAS, cgi);
}
if (cwi != GIMode::NONE) {
this->groupState().injection_control(group, Phase::WATER, cwi);
}
}
}
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
anyMSWellOpenLocal() const
{
for (const auto& well : wells_ecl_) {
if (well.isMultiSegment()) {
return true;
}
}
return false;
}
template<typename TypeTag>
const Well&
BlackoilWellModel<TypeTag>::
getWellEcl(const std::string& well_name) const
{
// finding the iterator of the well in wells_ecl
auto well_ecl = std::find_if(wells_ecl_.begin(),
wells_ecl_.end(),
[&well_name](const Well& elem)->bool {
return elem.name() == well_name;
});
assert(well_ecl != wells_ecl_.end());
return *well_ecl;
}
template<typename TypeTag>
typename BlackoilWellModel<TypeTag>::WellInterfacePtr
BlackoilWellModel<TypeTag>::
getWell(const std::string& well_name) const
{
// finding the iterator of the well in wells_ecl
auto well = std::find_if(well_container_.begin(),
well_container_.end(),
[&well_name](const WellInterfacePtr& elem)->bool {
return elem->name() == well_name;
});
assert(well != well_container_.end());
return *well;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateGroupIndividualControls(Opm::DeferredLogger& deferred_logger, std::set<std::string>& switched_groups)
{
const int reportStepIdx = ebosSimulator_.episodeIndex();
const int nupcol = schedule()[reportStepIdx].nupcol();
const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
// don't switch group control when iterationIdx > nupcol
// to avoid oscilations between group controls
if (iterationIdx > nupcol)
return;
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
updateGroupIndividualControl(fieldGroup, deferred_logger, switched_groups);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateGroupIndividualControl(const Group& group, Opm::DeferredLogger& deferred_logger, std::set<std::string>& switched_groups) {
const int reportStepIdx = ebosSimulator_.episodeIndex();
const bool skip = switched_groups.count(group.name());
if (!skip && group.isInjectionGroup())
{
const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS};
for (Phase phase : all) {
if (!group.hasInjectionControl(phase)) {
continue;
}
Group::InjectionCMode newControl = checkGroupInjectionConstraints(group, phase);
if (newControl != Group::InjectionCMode::NONE)
{
switched_groups.insert(group.name());
actionOnBrokenConstraints(group, newControl, phase, deferred_logger);
}
}
}
if (!skip && group.isProductionGroup()) {
Group::ProductionCMode newControl = checkGroupProductionConstraints(group, deferred_logger);
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
const auto controls = group.productionControls(summaryState);
if (newControl != Group::ProductionCMode::NONE)
{
switched_groups.insert(group.name());
actionOnBrokenConstraints(group, controls.exceed_action, newControl, deferred_logger);
}
}
// call recursively down the group hiearchy
for (const std::string& groupName : group.groups()) {
updateGroupIndividualControl( schedule().getGroup(groupName, reportStepIdx), deferred_logger, switched_groups);
}
}
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
checkGroupConstraints(const Group& group, Opm::DeferredLogger& deferred_logger) const {
const int reportStepIdx = ebosSimulator_.episodeIndex();
if (group.isInjectionGroup()) {
const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS};
for (Phase phase : all) {
if (!group.hasInjectionControl(phase)) {
continue;
}
Group::InjectionCMode newControl = checkGroupInjectionConstraints(group, phase);
if (newControl != Group::InjectionCMode::NONE) {
return true;
}
}
}
if (group.isProductionGroup()) {
Group::ProductionCMode newControl = checkGroupProductionConstraints(group, deferred_logger);
if (newControl != Group::ProductionCMode::NONE)
{
return true;
}
}
// call recursively down the group hiearchy
bool violated = false;
for (const std::string& groupName : group.groups()) {
violated = violated || checkGroupConstraints( schedule().getGroup(groupName, reportStepIdx), deferred_logger);
}
return violated;
}
template<typename TypeTag>
Group::ProductionCMode
BlackoilWellModel<TypeTag>::
checkGroupProductionConstraints(const Group& group, Opm::DeferredLogger& deferred_logger) const {
const int reportStepIdx = ebosSimulator_.episodeIndex();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
const auto& comm = ebosSimulator_.vanguard().grid().comm();
const auto& well_state = this->wellState();
const auto controls = group.productionControls(summaryState);
const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name());
if (group.has_control(Group::ProductionCMode::ORAT))
{
if (currentControl != Group::ProductionCMode::ORAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.oil_target < current_rate ) {
return Group::ProductionCMode::ORAT;
}
}
}
if (group.has_control(Group::ProductionCMode::WRAT))
{
if (currentControl != Group::ProductionCMode::WRAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.water_target < current_rate ) {
return Group::ProductionCMode::WRAT;
}
}
}
if (group.has_control(Group::ProductionCMode::GRAT))
{
if (currentControl != Group::ProductionCMode::GRAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.gas_target < current_rate ) {
return Group::ProductionCMode::GRAT;
}
}
}
if (group.has_control(Group::ProductionCMode::LRAT))
{
if (currentControl != Group::ProductionCMode::LRAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.liquid_target < current_rate ) {
return Group::ProductionCMode::LRAT;
}
}
}
if (group.has_control(Group::ProductionCMode::CRAT))
{
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "CRAT control for production groups not implemented" , deferred_logger);
}
if (group.has_control(Group::ProductionCMode::RESV))
{
if (currentControl != Group::ProductionCMode::RESV)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true);
current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true);
current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.resv_target < current_rate ) {
return Group::ProductionCMode::RESV;
}
}
}
if (group.has_control(Group::ProductionCMode::PRBL))
{
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "PRBL control for production groups not implemented", deferred_logger);
}
return Group::ProductionCMode::NONE;
}
template<typename TypeTag>
Group::InjectionCMode
BlackoilWellModel<TypeTag>::
checkGroupInjectionConstraints(const Group& group, const Phase& phase) const {
const int reportStepIdx = ebosSimulator_.episodeIndex();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
const auto& comm = ebosSimulator_.vanguard().grid().comm();
const auto& well_state = this->wellState();
int phasePos;
if (phase == Phase::GAS && phase_usage_.phase_used[BlackoilPhases::Vapour] )
phasePos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
else if (phase == Phase::OIL && phase_usage_.phase_used[BlackoilPhases::Liquid])
phasePos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
else if (phase == Phase::WATER && phase_usage_.phase_used[BlackoilPhases::Aqua] )
phasePos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
else
OPM_THROW(std::runtime_error, "Unknown phase" );
const auto& controls = group.injectionControls(phase, summaryState);
auto currentControl = this->groupState().injection_control(group.name(), phase);
if (controls.has_control(Group::InjectionCMode::RATE))
{
if (currentControl != Group::InjectionCMode::RATE)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.surface_max_rate < current_rate) {
return Group::InjectionCMode::RATE;
}
}
}
if (controls.has_control(Group::InjectionCMode::RESV))
{
if (currentControl != Group::InjectionCMode::RESV)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.resv_max_rate < current_rate) {
return Group::InjectionCMode::RESV;
}
}
}
if (controls.has_control(Group::InjectionCMode::REIN))
{
if (currentControl != Group::InjectionCMode::REIN)
{
double production_Rate = 0.0;
const Group& groupRein = schedule().getGroup(controls.reinj_group, reportStepIdx);
production_Rate += WellGroupHelpers::sumWellRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false);
// sum over all nodes
production_Rate = comm.sum(production_Rate);
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
// sum over all nodes
current_rate = comm.sum(current_rate);
if (controls.target_reinj_fraction*production_Rate < current_rate) {
return Group::InjectionCMode::REIN;
}
}
}
if (controls.has_control(Group::InjectionCMode::VREP))
{
if (currentControl != Group::InjectionCMode::VREP)
{
double voidage_rate = 0.0;
const Group& groupVoidage = schedule().getGroup(controls.voidage_group, reportStepIdx);
voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
// sum over all nodes
voidage_rate = comm.sum(voidage_rate);
double total_rate = 0.0;
total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true);
total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true);
total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true);
// sum over all nodes
total_rate = comm.sum(total_rate);
if (controls.target_void_fraction*voidage_rate < total_rate) {
return Group::InjectionCMode::VREP;
}
}
}
return Group::InjectionCMode::NONE;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
checkGconsaleLimits(const Group& group, WellState& well_state, Opm::DeferredLogger& deferred_logger)
{
const int reportStepIdx = ebosSimulator_.episodeIndex();
// call recursively down the group hiearchy
for (const std::string& groupName : group.groups()) {
checkGconsaleLimits( schedule().getGroup(groupName, reportStepIdx), well_state, deferred_logger);
}
// only for groups with gas injection controls
if (!group.hasInjectionControl(Phase::GAS)) {
return;
}
// check if gconsale is used for this group
if (!schedule()[reportStepIdx].gconsale().has(group.name()))
return;
std::ostringstream ss;
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
const auto& comm = ebosSimulator_.vanguard().grid().comm();
const auto& gconsale = schedule()[reportStepIdx].gconsale().get(group.name(), summaryState);
const Group::ProductionCMode& oldProductionControl = this->groupState().production_control(group.name());
int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
double production_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/false);
double injection_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/true);
// sum over all nodes
injection_rate = comm.sum(injection_rate);
production_rate = comm.sum(production_rate);
double sales_rate = production_rate - injection_rate;
double production_target = gconsale.sales_target + injection_rate;
// add import rate and substract consumption rate for group for gas
if (schedule()[reportStepIdx].gconsump().has(group.name())) {
const auto& gconsump = schedule()[reportStepIdx].gconsump().get(group.name(), summaryState);
if (phase_usage_.phase_used[BlackoilPhases::Vapour]) {
sales_rate += gconsump.import_rate;
sales_rate -= gconsump.consumption_rate;
production_target -= gconsump.import_rate;
production_target += gconsump.consumption_rate;
}
}
if (sales_rate > gconsale.max_sales_rate) {
switch(gconsale.max_proc) {
case GConSale::MaxProcedure::NONE: {
if (oldProductionControl != Group::ProductionCMode::GRAT && oldProductionControl != Group::ProductionCMode::NONE) {
ss << "Group sales exceed maximum limit, but the action is NONE for " + group.name() + ". Nothing happens";
}
break;
}
case GConSale::MaxProcedure::CON: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit CON not implemented", deferred_logger);
break;
}
case GConSale::MaxProcedure::CON_P: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit CON_P not implemented", deferred_logger);
break;
}
case GConSale::MaxProcedure::WELL: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit WELL not implemented", deferred_logger);
break;
}
case GConSale::MaxProcedure::PLUG: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit PLUG not implemented", deferred_logger);
break;
}
case GConSale::MaxProcedure::MAXR: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit MAXR not implemented", deferred_logger);
break;
}
case GConSale::MaxProcedure::END: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit END not implemented", deferred_logger);
break;
}
case GConSale::MaxProcedure::RATE: {
this->groupState().production_control(group.name(), Group::ProductionCMode::GRAT);
ss << "Maximum GCONSALE limit violated for " << group.name() << ". The group is switched from ";
ss << Group::ProductionCMode2String(oldProductionControl) << " to " << Group::ProductionCMode2String(Group::ProductionCMode::GRAT);
ss << " and limited by the maximum sales rate after consumption and import are considered" ;
this->groupState().update_grat_sales_target(group.name(), production_target);
break;
}
default:
throw("Invalid procedure for maximum rate limit selected for group" + group.name());
}
}
if (sales_rate < gconsale.min_sales_rate) {
const Group::ProductionCMode& currentProductionControl = this->groupState().production_control(group.name());
if ( currentProductionControl == Group::ProductionCMode::GRAT ) {
ss << "Group " + group.name() + " has sale rate less then minimum permitted value and is under GRAT control. \n";
ss << "The GRAT is increased to meet the sales minimum rate. \n";
this->groupState().update_grat_sales_target(group.name(), production_target);
//} else if () {//TODO add action for WGASPROD
//} else if () {//TODO add action for drilling queue
} else {
ss << "Group " + group.name() + " has sale rate less then minimum permitted value but cannot increase the group production rate \n";
ss << "or adjust gas production using WGASPROD or drill new wells to meet the sales target. \n";
ss << "Note that WGASPROD and drilling queues are not implemented in Flow. No action is taken. \n ";
}
}
if (gconsale.sales_target < 0.0) {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + " has sale rate target less then zero. Not implemented in Flow" , deferred_logger);
}
auto cc = Dune::MPIHelper::getCollectiveCommunication();
if (!ss.str().empty() && cc.rank() == 0)
deferred_logger.info(ss.str());
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
actionOnBrokenConstraints(const Group& group, const Group::ExceedAction& exceed_action, const Group::ProductionCMode& newControl, Opm::DeferredLogger& deferred_logger) {
const Group::ProductionCMode oldControl = this->groupState().production_control(group.name());
std::ostringstream ss;
switch(exceed_action) {
case Group::ExceedAction::NONE: {
if (oldControl != newControl && oldControl != Group::ProductionCMode::NONE) {
ss << "Group production exceed action is NONE for group " + group.name() + ". Nothing happens.";
}
break;
}
case Group::ExceedAction::CON: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit CON not implemented", deferred_logger);
break;
}
case Group::ExceedAction::CON_PLUS: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit CON_PLUS not implemented", deferred_logger);
break;
}
case Group::ExceedAction::WELL: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit WELL not implemented", deferred_logger);
break;
}
case Group::ExceedAction::PLUG: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit PLUG not implemented", deferred_logger);
break;
}
case Group::ExceedAction::RATE: {
if (oldControl != newControl) {
this->groupState().production_control(group.name(), newControl);
ss << "Switching production control mode for group "<< group.name()
<< " from " << Group::ProductionCMode2String(oldControl)
<< " to " << Group::ProductionCMode2String(newControl);
}
break;
}
default:
throw("Invalid procedure for maximum rate limit selected for group" + group.name());
}
auto cc = Dune::MPIHelper::getCollectiveCommunication();
if (!ss.str().empty() && cc.rank() == 0)
deferred_logger.info(ss.str());
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
actionOnBrokenConstraints(const Group& group, const Group::InjectionCMode& newControl, const Phase& controlPhase, Opm::DeferredLogger& deferred_logger) {
auto oldControl = this->groupState().injection_control(group.name(), controlPhase);
std::ostringstream ss;
if (oldControl != newControl) {
const std::string from = Group::InjectionCMode2String(oldControl);
ss << "Switching injection control mode for group "<< group.name()
<< " from " << Group::InjectionCMode2String(oldControl)
<< " to " << Group::InjectionCMode2String(newControl);
this->groupState().injection_control(group.name(), controlPhase, newControl);
}
auto cc = Dune::MPIHelper::getCollectiveCommunication();
if (!ss.str().empty() && cc.rank() == 0)
deferred_logger.info(ss.str());
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateGroupHigherControls(Opm::DeferredLogger& deferred_logger, std::set<std::string>& switched_groups)
{
const int reportStepIdx = ebosSimulator_.episodeIndex();
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
checkGroupHigherConstraints(fieldGroup, deferred_logger, switched_groups);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
checkGroupHigherConstraints(const Group& group, Opm::DeferredLogger& deferred_logger, std::set<std::string>& switched_groups)
{
// Set up coefficients for RESV <-> surface rate conversion.
// Use the pvtRegionIdx from the top cell of the first well.
// TODO fix this!
// This is only used for converting RESV rates.
// What is the proper approach?
const auto& comm = ebosSimulator_.vanguard().grid().comm();
const int fipnum = 0;
int pvtreg = well_perf_data_.empty() || well_perf_data_[0].empty()
? pvt_region_idx_[0]
: pvt_region_idx_[well_perf_data_[0][0].cell_index];
if ( comm.size() > 1)
{
// Just like in the sequential case the pvtregion is determined
// by the first cell of the first well. What is the first well
// is decided by the order in the Schedule using Well::seqIndex()
int firstWellIndex = well_perf_data_.empty() ?
std::numeric_limits<int>::max() : wells_ecl_[0].seqIndex();
auto regIndexPair = std::make_pair(pvtreg, firstWellIndex);
std::vector<decltype(regIndexPair)> pairs(comm.size());
comm.allgather(&regIndexPair, 1, pairs.data());
pvtreg = std::min_element(pairs.begin(), pairs.end(),
[](const auto& p1, const auto& p2){ return p1.second < p2.second;})
->first;
}
std::vector<double> resv_coeff(phase_usage_.num_phases, 0.0);
rateConverter_->calcCoeff(fipnum, pvtreg, resv_coeff);
const int reportStepIdx = ebosSimulator_.episodeIndex();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
std::vector<double> rates(phase_usage_.num_phases, 0.0);
const bool skip = switched_groups.count(group.name()) || group.name() == "FIELD";
if (!skip && group.isInjectionGroup()) {
// Obtain rates for group.
for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) {
const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ true);
// Sum over all processes
rates[phasePos] = comm.sum(local_current_rate);
}
const Phase all[] = { Phase::WATER, Phase::OIL, Phase::GAS };
for (Phase phase : all) {
// Check higher up only if under individual (not FLD) control.
auto currentControl = this->groupState().injection_control(group.name(), phase);
if (currentControl != Group::InjectionCMode::FLD && group.injectionGroupControlAvailable(phase)) {
const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx);
const std::pair<bool, double> changed = WellGroupHelpers::checkGroupConstraintsInj(
group.name(),
group.parent(),
parentGroup,
this->wellState(),
this->groupState(),
reportStepIdx,
guideRate_.get(),
rates.data(),
phase,
phase_usage_,
group.getGroupEfficiencyFactor(),
schedule(),
summaryState,
resv_coeff,
deferred_logger);
if (changed.first) {
switched_groups.insert(group.name());
actionOnBrokenConstraints(group, Group::InjectionCMode::FLD, phase, deferred_logger);
}
}
}
}
if (!skip && group.isProductionGroup()) {
// Obtain rates for group.
for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) {
const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ false);
// Sum over all processes
rates[phasePos] = -comm.sum(local_current_rate);
}
// Check higher up only if under individual (not FLD) control.
const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name());
if (currentControl != Group::ProductionCMode::FLD && group.productionGroupControlAvailable()) {
const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx);
const std::pair<bool, double> changed = WellGroupHelpers::checkGroupConstraintsProd(
group.name(),
group.parent(),
parentGroup,
this->wellState(),
this->groupState(),
reportStepIdx,
guideRate_.get(),
rates.data(),
phase_usage_,
group.getGroupEfficiencyFactor(),
schedule(),
summaryState,
resv_coeff,
deferred_logger);
if (changed.first) {
switched_groups.insert(group.name());
const auto exceed_action = group.productionControls(summaryState).exceed_action;
actionOnBrokenConstraints(group, exceed_action, Group::ProductionCMode::FLD, deferred_logger);
}
}
}
// call recursively down the group hiearchy
for (const std::string& groupName : group.groups()) {
checkGroupHigherConstraints( schedule().getGroup(groupName, reportStepIdx), deferred_logger, switched_groups);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateEclWells(const int timeStepIdx, const std::unordered_set<std::string>& wells) {
const auto& schedule = this->ebosSimulator_.vanguard().schedule();
for (const auto& wname : wells) {
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()) {
auto well_index = std::distance( this->wells_ecl_.begin(), well_iter );
this->wells_ecl_[well_index] = schedule.getWell(wname, timeStepIdx);
const auto& well = this->wells_ecl_[well_index];
auto& pd = this->well_perf_data_[well_index];
auto pdIter = pd.begin();
for (const auto& conn : well.getConnections()) {
if (conn.state() != Connection::State::SHUT) {
pdIter->connection_transmissibility_factor = conn.CF();
++pdIter;
}
}
this->wellState().updateStatus(well_index, well.getStatus());
this->wellState().resetConnectionTransFactors(well_index, pd);
this->prod_index_calc_[well_index].reInit(well);
}
}
}
template<typename TypeTag>
double
BlackoilWellModel<TypeTag>::
wellPI(const int well_index) const
{
const auto& pu = this->phase_usage_;
const auto np = this->numPhases();
const auto* pi = &this->wellState().productivityIndex()[np*well_index + 0];
const auto preferred = this->wells_ecl_[well_index].getPreferredPhase();
switch (preferred) { // Should really have LIQUID = OIL + WATER here too...
case Phase::WATER:
return pu.phase_used[BlackoilPhases::PhaseIndex::Aqua]
? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Aqua]]
: 0.0;
case Phase::OIL:
return pu.phase_used[BlackoilPhases::PhaseIndex::Liquid]
? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Liquid]]
: 0.0;
case Phase::GAS:
return pu.phase_used[BlackoilPhases::PhaseIndex::Vapour]
? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Vapour]]
: 0.0;
default:
throw std::invalid_argument {
"Unsupported preferred phase " +
std::to_string(static_cast<int>(preferred))
};
}
}
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_index = std::distance(this->wells_ecl_.begin(), well_iter);
return this->wellPI(well_index);
}
template <typename TypeTag>
int
BlackoilWellModel<TypeTag>::
reportStepIndex() const
{
return std::max(this->ebosSimulator_.episodeIndex(), 0);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger)
{
if (this->last_run_wellpi_.has_value() && (*this->last_run_wellpi_ == timeStepIdx)) {
// We've already run WELPI scaling for this report step. Most
// common for the very first report step. Don't redo WELPI scaling.
return;
}
auto hasWellPIEvent = [this, timeStepIdx](const int well_index) -> bool
{
return this->schedule()[timeStepIdx].wellgroup_events()
.hasEvent(this->wells_ecl_[well_index].name(),
ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX);
};
auto updateEclWell = [this, timeStepIdx](const int well_index) -> void
{
const auto& schedule = this->schedule();
const auto& wname = this->wells_ecl_[well_index].name();
this->wells_ecl_[well_index] = schedule.getWell(wname, timeStepIdx);
const auto& well = this->wells_ecl_[well_index];
auto& pd = this->well_perf_data_[well_index];
auto pdIter = pd.begin();
for (const auto& conn : well.getConnections()) {
if (conn.state() != Connection::State::SHUT) {
pdIter->connection_transmissibility_factor = conn.CF();
++pdIter;
}
}
this->wellState().resetConnectionTransFactors(well_index, pd);
this->prod_index_calc_[well_index].reInit(well);
};
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);
};
// Minimal well setup to compute PI/II values
{
auto saved_previous_wgstate = this->prevWGState();
this->commitWGState();
well_container_ = createWellContainer(timeStepIdx);
for (auto& well : well_container_) {
well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_);
}
std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
for (auto& well : well_container_) {
well->updatePerforatedCell(is_cell_perforated_);
}
this->calculateProductivityIndexValues(local_deferredLogger);
this->commitWGState(std::move(saved_previous_wgstate));
}
const auto nw = this->numLocalWells();
for (auto wellID = 0*nw; wellID < nw; ++wellID) {
if (hasWellPIEvent(wellID)) {
rescaleWellPI(wellID, this->wellPI(wellID));
updateEclWell(wellID);
}
}
this->last_run_wellpi_ = timeStepIdx;
}
template <typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
wasDynamicallyShutThisTimeStep(const int well_index) const
{
return this->closed_this_step_.find(this->wells_ecl_[well_index].name())
!= this->closed_this_step_.end();
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, const WellStateFullyImplicitBlackoil& wellState) {
for (const std::string& groupName : group.groups()) {
const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
updateWsolvent(groupTmp, schedule, reportStepIdx, wellState);
}
if (group.isProductionGroup())
return;
auto currentGroupControl = this->groupState().injection_control(group.name(), Phase::GAS);
if( currentGroupControl == Group::InjectionCMode::REIN ) {
int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
const auto& controls = group.injectionControls(Phase::GAS, summaryState);
const Group& groupRein = schedule.getGroup(controls.reinj_group, reportStepIdx);
double gasProductionRate = WellGroupHelpers::sumWellRates(groupRein, schedule, wellState, reportStepIdx, gasPos, /*isInjector*/false);
double solventProductionRate = WellGroupHelpers::sumSolventRates(groupRein, schedule, wellState, reportStepIdx, /*isInjector*/false);
const auto& comm = ebosSimulator_.vanguard().grid().comm();
solventProductionRate = comm.sum(solventProductionRate);
gasProductionRate = comm.sum(gasProductionRate);
double wsolvent = 0.0;
if (std::abs(gasProductionRate) > 1e-6)
wsolvent = solventProductionRate / gasProductionRate;
setWsolvent(group, schedule, reportStepIdx, wsolvent);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent) {
for (const std::string& groupName : group.groups()) {
const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
setWsolvent(groupTmp, schedule, reportStepIdx, wsolvent);
}
for (const std::string& wellName : group.wells()) {
const auto& wellTmp = schedule.getWell(wellName, reportStepIdx);
if (wellTmp.getStatus() == Well::Status::SHUT)
continue;
auto well = getWell(wellName);
well->setWsolvent(wsolvent);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignWellGuideRates(data::Wells& wsrpt) const
{
for (const auto& well : this->wells_ecl_) {
auto xwPos = wsrpt.find(well.name());
if (xwPos == wsrpt.end()) { // No well results. Unexpected.
continue;
}
xwPos->second.guide_rates = this->getGuideRateValues(well);
}
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignShutConnections(data::Wells& wsrpt) const
{
auto wellID = 0;
for (const auto& well : this->wells_ecl_) {
auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<>
xwel.dynamicStatus = this->schedule()
.getWell(well.name(), this->reportStepIndex()).getStatus();
const auto wellIsOpen = xwel.dynamicStatus == Well::Status::OPEN;
auto skip = [wellIsOpen](const Connection& conn)
{
return wellIsOpen && (conn.state() != Connection::State::SHUT);
};
if (this->wellTestState_.hasWellClosed(well.name()) &&
!this->wasDynamicallyShutThisTimeStep(wellID))
{
xwel.dynamicStatus = well.getAutomaticShutIn()
? Well::Status::SHUT : Well::Status::STOP;
}
auto& xcon = xwel.connections;
for (const auto& conn : well.getConnections()) {
if (skip(conn)) {
continue;
}
auto& xc = xcon.emplace_back();
xc.index = conn.global_index();
xc.pressure = xc.reservoir_rate = 0.0;
xc.effective_Kh = conn.Kh();
xc.trans_factor = conn.CF();
}
++wellID;
}
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignGroupValues(const int reportStepIdx,
const Schedule& sched,
std::map<std::string, data::GroupData>& gvalues) const
{
const auto groupGuideRates =
this->calculateAllGroupGuiderates(reportStepIdx, sched);
for (const auto& gname : sched.groupNames(reportStepIdx)) {
const auto& grup = sched.getGroup(gname, reportStepIdx);
auto& gdata = gvalues[gname];
this->assignGroupControl(grup, gdata);
this->assignGroupGuideRates(grup, groupGuideRates, gdata);
}
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignNodeValues(std::map<std::string, data::NodeData>& nodevalues) const
{
nodevalues.clear();
for (const auto& [node, pressure] : node_pressures_) {
nodevalues.emplace(node, data::NodeData{pressure});
}
}
template<typename TypeTag>
std::unordered_map<std::string, data::GroupGuideRates>
BlackoilWellModel<TypeTag>::
calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const
{
auto gr = std::unordered_map<std::string, data::GroupGuideRates>{};
auto up = std::vector<std::string>{};
// Start at well level, accumulate contributions towards root of
// group tree (FIELD group).
for (const auto& wname : sched.wellNames(reportStepIdx)) {
if (! (this->wellState().hasWellRates(wname) &&
this->guideRate_->has(wname)))
{
continue;
}
const auto& well = sched.getWell(wname, reportStepIdx);
const auto& parent = well.groupName();
if (parent == "FIELD") {
// Well parented directly to "FIELD". Inadvisable and
// unexpected, but nothing to do about that here. Just skip
// this guide rate contribution.
continue;
}
auto& grval = well.isInjector()
? gr[parent].injection
: gr[parent].production;
grval += this->getGuideRateValues(well);
up.push_back(parent);
}
// Propagate accumulated guide rates up towards root of group tree.
// Override accumulation if there is a GUIDERAT specification that
// applies to a group.
std::sort(up.begin(), up.end());
auto start = 0*up.size();
auto u = std::unique(up.begin(), up.end());
auto nu = std::distance(up.begin(), u);
while (nu > 0) {
const auto ntot = up.size();
for (auto gi = 0*nu; gi < nu; ++gi) {
const auto& gname = up[start + gi];
const auto& group = sched.getGroup(gname, reportStepIdx);
if (this->guideRate_->has(gname)) {
gr[gname].production = this->getGuideRateValues(group);
}
if (this->guideRate_->has(gname, Opm::Phase::WATER)
|| this->guideRate_->has(gname, Opm::Phase::GAS)) {
gr[gname].injection = this->getGuideRateInjectionGroupValues(group);
}
const auto parent = group.parent();
if (parent == "FIELD") { continue; }
gr[parent].injection += gr[gname].injection;
gr[parent].production += gr[gname].production;
up.push_back(parent);
}
start = ntot;
auto begin = up.begin() + ntot;
std::sort(begin, up.end());
u = std::unique(begin, up.end());
nu = std::distance(begin, u);
}
return gr;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignGroupControl(const Group& group, data::GroupData& gdata) const
{
const auto& gname = group.name();
const auto grup_type = group.getGroupType();
auto& cgc = gdata.currentControl;
cgc.currentProdConstraint =
::Opm::Group::ProductionCMode::NONE;
cgc.currentGasInjectionConstraint =
cgc.currentWaterInjectionConstraint =
::Opm::Group::InjectionCMode::NONE;
if (this->groupState().has_production_control(gname)) {
cgc.currentProdConstraint = this->groupState().production_control(gname);
}
if ((grup_type == ::Opm::Group::GroupType::INJECTION) ||
(grup_type == ::Opm::Group::GroupType::MIXED))
{
if (this->groupState().has_injection_control(gname, ::Opm::Phase::WATER)) {
cgc.currentWaterInjectionConstraint = this->groupState().injection_control(gname, Phase::WATER);
}
if (this->groupState().has_injection_control(gname, ::Opm::Phase::GAS)) {
cgc.currentGasInjectionConstraint = this->groupState().injection_control(gname, Phase::GAS);
}
}
}
template <typename TypeTag>
data::GuideRateValue
BlackoilWellModel<TypeTag>::
getGuideRateValues(const Well& well) const
{
auto grval = data::GuideRateValue{};
assert (this->guideRate_ != nullptr);
const auto& wname = well.name();
if (! this->wellState().hasWellRates(wname)) {
// No flow rates for 'wname' -- might be before well comes
// online (e.g., for the initial condition before simulation
// starts).
return grval;
}
if (! this->guideRate_->has(wname)) {
// No guiderates exist for 'wname'.
return grval;
}
const auto qs = WellGroupHelpers::
getWellRateVector(this->wellState(), this->phase_usage_, wname);
this->getGuideRateValues(qs, well.isInjector(), wname, grval);
return grval;
}
template <typename TypeTag>
data::GuideRateValue
BlackoilWellModel<TypeTag>::
getGuideRateInjectionGroupValues(const Group& group) const
{
auto grval = data::GuideRateValue{};
assert (this->guideRate_ != nullptr);
const auto& gname = group.name();
if (this->guideRate_->has(gname, Opm::Phase::GAS)) {
grval.set(data::GuideRateValue::Item::Gas,
this->guideRate_->get(gname, Opm::Phase::GAS));
}
if (this->guideRate_->has(gname, Opm::Phase::WATER)) {
grval.set(data::GuideRateValue::Item::Water,
this->guideRate_->get(gname, Opm::Phase::WATER));
}
return grval;
}
template <typename TypeTag>
data::GuideRateValue
BlackoilWellModel<TypeTag>::
getGuideRateValues(const Group& group) const
{
auto grval = data::GuideRateValue{};
assert (this->guideRate_ != nullptr);
const auto& gname = group.name();
if ( ! this->groupState().has_production_rates(gname)) {
// No flow rates for production group 'gname' -- might be before group comes
// online (e.g., for the initial condition before simulation
// starts).
return grval;
}
if (! this->guideRate_->has(gname)) {
// No guiderates exist for 'gname'.
return grval;
}
const auto qs = WellGroupHelpers::getProductionGroupRateVector(this->groupState(), this->phase_usage_, gname);
const auto is_inj = false; // This procedure only applies to G*PGR.
this->getGuideRateValues(qs, is_inj, gname, grval);
return grval;
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
getGuideRateValues(const GuideRate::RateVector& qs,
const bool is_inj,
const std::string& wgname,
data::GuideRateValue& grval) const
{
auto getGR = [this, &wgname, &qs](const GuideRateModel::Target t)
{
return this->guideRate_->get(wgname, t, qs);
};
// Note: GuideRate does currently (2020-07-20) not support Target::RES.
grval.set(data::GuideRateValue::Item::Gas,
getGR(GuideRateModel::Target::GAS));
grval.set(data::GuideRateValue::Item::Water,
getGR(GuideRateModel::Target::WAT));
if (! is_inj) {
// Producer. Extract "all" guiderate values.
grval.set(data::GuideRateValue::Item::Oil,
getGR(GuideRateModel::Target::OIL));
}
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assignGroupGuideRates(const Group& group,
const std::unordered_map<std::string, data::GroupGuideRates>& groupGuideRates,
data::GroupData& gdata) const
{
auto& prod = gdata.guideRates.production; prod.clear();
auto& inj = gdata.guideRates.injection; inj .clear();
auto xgrPos = groupGuideRates.find(group.name());
if ((xgrPos == groupGuideRates.end()) ||
!this->guideRate_->has(group.name()))
{
// No guiderates defined for this group.
return;
}
const auto& xgr = xgrPos->second;
prod = xgr.production;
inj = xgr.injection;
}
template <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
computeWellTemperature()
{
if (!has_energy_)
return;
int np = numPhases();
double cellInternalEnergy;
double cellBinv;
double cellDensity;
double perfPhaseRate;
const int nw = numLocalWells();
for (auto wellID = 0*nw; wellID < nw; ++wellID) {
const Well& well = wells_ecl_[wellID];
if (well.isInjector())
continue;
int connpos = 0;
for (int i = 0; i < wellID; ++i) {
connpos += well_perf_data_[i].size();
}
connpos *= np;
double weighted_temperature = 0.0;
double total_weight = 0.0;
auto& well_info = *local_parallel_well_info_[wellID];
const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size());
for (int perf = 0; perf < num_perf_this_well; ++perf) {
const int cell_idx = well_perf_data_[wellID][perf].cell_index;
const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& fs = intQuants.fluidState();
double cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
double weight_factor = 0.0;
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
cellInternalEnergy = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
cellBinv = fs.invB(phaseIdx).value();
cellDensity = fs.density(phaseIdx).value();
perfPhaseRate = this->wellState().perfPhaseRates()[connpos + perf*np + phaseIdx ];
weight_factor += cellDensity * perfPhaseRate/cellBinv * cellInternalEnergy/cellTemperatures;
}
total_weight += weight_factor;
weighted_temperature += weight_factor * cellTemperatures;
}
weighted_temperature = well_info.communication().sum(weighted_temperature);
total_weight = well_info.communication().sum(total_weight);
this->wellState().temperature()[wellID] = weighted_temperature/total_weight;
}
}
} // namespace Opm