opm-simulators/opm/simulators/wells/BlackoilWellModel_impl.hpp
Atgeirr Flø Rasmussen bddeaba880 Delete obsolete include directives and functions.
Also make other minor adjustments, in particular moving the function
getCubeDim() from WellsManagerDetail to WellHelpers.hpp.
2019-11-25 10:58:44 +01:00

1931 lines
76 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/simulators/wells/SimFIBODetails.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
namespace Opm {
template<typename TypeTag>
BlackoilWellModel<TypeTag>::
BlackoilWellModel(Simulator& ebosSimulator)
: ebosSimulator_(ebosSimulator)
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
{
terminal_output_ = false;
if (ebosSimulator.gridView().comm().rank() == 0)
terminal_output_ = EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput);
// Create the guide rate container.
guideRate_.reset(new GuideRate (ebosSimulator_.vanguard().schedule()));
// calculate the number of elements of the compressed sequential grid. this needs
// to be done in two steps because the dune communicator expects a reference as
// argument for sum()
const auto& gridView = ebosSimulator_.gridView();
number_of_cells_ = gridView.size(/*codim=*/0);
global_nc_ = gridView.comm().sum(number_of_cells_);
// Set up cartesian mapping.
const auto& grid = ebosSimulator_.vanguard().grid();
const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid),
cartDims[0]*cartDims[1]*cartDims[2]);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
init()
{
const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState();
extractLegacyCellPvtRegionIndex_();
extractLegacyDepth_();
phase_usage_ = phaseUsageFromDeck(eclState);
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(number_of_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();
const auto& cartesianSize = Opm::UgGridHelpers::cartDims(grid());
// 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 i = connection.getI();
int j = connection.getJ();
int k = connection.getK();
int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
int compressed_idx = cartesian_to_compressed_.at(cart_grid_idx);
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>
void
BlackoilWellModel<TypeTag>::
beginReportStep(const int timeStepIdx)
{
Opm::DeferredLogger local_deferredLogger;
const Grid& grid = ebosSimulator_.vanguard().grid();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
// Make wells_ecl_ contain only this partition's non-shut wells.
{
const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
auto is_shut_or_defunct = [&defunct_well_names](const Well& well) {
return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end());
};
auto w = schedule().getWells(timeStepIdx);
w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end());
wells_ecl_.swap(w);
}
initializeWellPerfData();
// 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_);
// The well state initialize bhp with the cell pressure in the top cell.
// We must therefore provide it with updated cell pressures
size_t nc = number_of_cells_;
std::vector<double> cellPressures(nc, 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)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity) {
continue;
}
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// copy of get perfpressure in Standard well
// exept for value
double perf_pressure = 0.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();
}
}
cellPressures[cellIdx] = perf_pressure;
}
well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_, well_perf_data_, summaryState);
// handling MS well related
if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
well_state_.initWellStateMSWell(wells_ecl_, phase_usage_, &previous_well_state_);
}
const int nw = wells_ecl_.size();
for (int w = 0; w <nw; ++w) {
const auto& well = wells_ecl_[w];
const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+ ScheduleEvents::PRODUCTION_UPDATE
+ ScheduleEvents::INJECTION_UPDATE
+ ScheduleEvents::NEW_WELL;
if(!schedule().hasWellEvent(well.name(), effective_events_mask, timeStepIdx))
continue;
if (well.isProducer()) {
const auto controls = well.productionControls(summaryState);
well_state_.currentProductionControls()[w] = controls.cmode;
}
else {
const auto controls = well.injectionControls(summaryState);
well_state_.currentInjectionControls()[w] = controls.cmode;
}
}
const Group& fieldGroup = schedule().getGroup("FIELD", timeStepIdx);
wellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, well_state_);
// Compute reservoir volumes for RESV controls.
rateConverter_.reset(new RateConverterType (phase_usage_,
std::vector<int>(number_of_cells_, 0)));
rateConverter_->template defineState<ElementContext>(ebosSimulator_);
// update VFP properties
vfp_properties_.reset (new VFPProperties<VFPInjProperties,VFPProdProperties> (
schedule().getVFPInjTables(timeStepIdx),
schedule().getVFPProdTables(timeStepIdx)) );
// update the previous well state. This is used to restart failed steps.
previous_well_state_ = well_state_;
}
// called at the beginning of a time step
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
beginTimeStep() {
Opm::DeferredLogger local_deferredLogger;
well_state_ = previous_well_state_;
const int reportStepIdx = ebosSimulator_.episodeIndex();
const double simulationTime = ebosSimulator_.time();
int exception_thrown = 0;
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_, number_of_cells_);
}
// 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 (has_polymer_)
{
const Grid& grid = ebosSimulator_.vanguard().grid();
if (PolymerModule::hasPlyshlog() || GET_PROP_VALUE(TypeTag, EnablePolymerMW) ) {
computeRepRadiusPerfLength(grid, local_deferredLogger);
}
}
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(local_deferredLogger, exception_thrown, "beginTimeStep() failed.", 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);
}
//compute well guideRates
const int np = numPhases();
for (const auto& well : well_container_) {
const auto wpot = well_state_.wellPotentials().data() + well->indexOfWell()*np;
const double oilpot = (phase_usage_.phase_used[BlackoilPhases::Liquid] > 0)
? wpot[phase_usage_.phase_pos[BlackoilPhases::Liquid]]
: 0.0;
const double gaspot = (phase_usage_.phase_used[BlackoilPhases::Vapour] > 0)
? wpot[phase_usage_.phase_pos[BlackoilPhases::Vapour]]
: 0.0;
const double waterpot = (phase_usage_.phase_used[BlackoilPhases::Aqua] > 0)
? wpot[phase_usage_.phase_pos[BlackoilPhases::Aqua]]
: 0.0;
guideRate_->compute(well->name(), reportStepIdx, simulationTime, oilpot, gaspot, waterpot);
}
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
wellGroupHelpers::updateGuideRateForGroups(fieldGroup, schedule(), phase_usage_, reportStepIdx, simulationTime, guideRate_.get(), well_state_);
// compute wsolvent fraction for REIN wells
updateWsolvent(fieldGroup, schedule(), reportStepIdx, well_state_);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx, const double simulationTime, Opm::DeferredLogger& deferred_logger) {
const auto& wtest_config = schedule().wtestConfig(timeStepIdx);
if (wtest_config.size() != 0) { // there is a WTEST request
// average B factors are required for the convergence checking of well equations
// Note: this must be done on all processes, even those with
// no wells needing testing, otherwise we will have locking.
std::vector< Scalar > B_avg(numComponents(), Scalar() );
computeAverageFormationFactor(B_avg);
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_, number_of_cells_);
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_, B_avg, simulationTime, timeStepIdx,
testing_reason, well_state_, wellTestState_, deferred_logger);
}
}
}
// called at the end of a report step
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
endReportStep() {
}
// called at the end of a report step
template<typename TypeTag>
const SimulatorReport&
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) {
Opm::DeferredLogger local_deferredLogger;
for (const auto& well : well_container_) {
if (GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well->isInjector()) {
well->updateWaterThroughput(dt, well_state_);
}
}
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;
const int reportStepIdx = ebosSimulator_.episodeIndex();
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);
}
previous_well_state_ = well_state_;
Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
if (terminal_output_) {
global_deferredLogger.logMessages();
}
}
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();
// Make wells_ecl_ contain only this partition's non-shut wells.
{
const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
auto is_shut_or_defunct = [&defunct_well_names](const Well& well) {
return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end());
};
auto w = schedule().getWells(report_step);
w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end());
wells_ecl_.swap(w);
}
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());
well_state_.resize(wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage, well_perf_data_, summaryState); // Resize for restart step
wellsToState(restartValues.wells, phaseUsage, handle_ms_well, well_state_);
}
// for ecl compatible restart the current controls are not written
const auto& ioCfg = eclState().getIOConfig();
const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST();
if (true || ecl_compatible_rst) { // always set the control from the schedule
for (int w = 0; w <nw; ++w) {
const auto& well = wells_ecl_[w];
if (well.isProducer()) {
const auto controls = well.productionControls(summaryState);
well_state_.currentProductionControls()[w] = controls.cmode;
}
else {
const auto controls = well.injectionControls(summaryState);
well_state_.currentInjectionControls()[w] = controls.cmode;
}
}
}
previous_well_state_ = well_state_;
initial_step_ = false;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initializeWellPerfData()
{
const auto& grid = ebosSimulator_.vanguard().grid();
const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
well_perf_data_.resize(wells_ecl_.size());
first_perf_index_.clear();
first_perf_index_.resize(wells_ecl_.size() + 1, 0);
int well_index = 0;
for (const auto& well : wells_ecl_) {
well_perf_data_[well_index].clear();
well_perf_data_[well_index].reserve(well.getConnections().size());
for (const auto& completion : well.getConnections()) {
if (completion.state() == Connection::State::OPEN) {
const int i = completion.getI();
const int j = completion.getJ();
const int k = completion.getK();
const int cart_grid_indx = i + cartDims[0] * (j + cartDims[1] * k);
const int active_index = cartesian_to_compressed_[cart_grid_indx];
if (active_index < 0) {
const std::string msg
= ("Cell with i,j,k indices " + std::to_string(i) + " " + std::to_string(j) + " "
+ std::to_string(k) + " not found in grid (well = " + well.name() + ").");
OPM_THROW(std::runtime_error, msg);
} else {
PerforationData pd;
pd.cell_index = active_index;
pd.connection_transmissibility_factor = completion.CF() * completion.wellPi();
pd.satnum_id = completion.satTableId();
well_perf_data_[well_index].push_back(pd);
}
} else {
if (completion.state() != Connection::State::SHUT) {
OPM_THROW(std::runtime_error,
"Completion state: " << Connection::State2String(completion.state()) << " not handled");
}
}
}
first_perf_index_[well_index + 1] = first_perf_index_[well_index] + well_perf_data_[well_index].size();
++well_index;
}
}
template<typename TypeTag>
std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
BlackoilWellModel<TypeTag>::
createWellContainer(const int time_step)
{
std::vector<WellInterfacePtr> well_container;
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();
// A new WCON keywords can re-open a well that was closed/shut due to Physical limit
if ( wellTestState_.hasWellClosed(well_name)) {
// TODO: more checking here, to make sure this standard more specific and complete
// maybe there is some WCON keywords will not open the well
if (well_state_.effectiveEventsOccurred(w)) {
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.
well_state_.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
// TODO: make a function from well_state side to handle the following
well_state_.thp()[w] = 0.;
well_state_.bhp()[w] = 0.;
const int np = numPhases();
for (int p = 0; p < np; ++p) {
well_state_.wellRates()[np * w + p] = 0.;
well_state_.wellReservoirRates()[np * w + p] = 0.;
}
continue;
} else {
// stopped wells are added to the container but marked as stopped
well_state_.thp()[w] = 0.;
wellIsStopped = true;
}
}
// Use the pvtRegionIdx from the top cell
const int well_cell_top = well_perf_data_[w][0].cell_index;
const int pvtreg = pvt_region_idx_[well_cell_top];
if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) {
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl,
time_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
w,
first_perf_index_[w],
well_perf_data_[w]));
} else {
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl,
time_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
w,
first_perf_index_[w],
well_perf_data_[w]));
}
if (wellIsStopped)
well_container.back()->stopWell();
}
}
return well_container;
}
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);
}
const Well& well_ecl = wells_ecl_[index_well_ecl];
// Use the pvtRegionIdx from the top cell
const int well_cell_top = well_perf_data_[index_well_ecl][0].cell_index;
const int pvtreg = pvt_region_idx_[well_cell_top];
if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) {
return WellInterfacePtr(new StandardWell<TypeTag>(well_ecl,
report_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
index_well_ecl,
first_perf_index_[index_well_ecl],
well_perf_data_[index_well_ecl]));
} else {
return WellInterfacePtr(new MultisegmentWell<TypeTag>(well_ecl,
report_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
index_well_ecl,
first_perf_index_[index_well_ecl],
well_perf_data_[index_well_ecl]));
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assemble(const int iterationIdx,
const double dt)
{
last_report_ = SimulatorReport();
if ( ! wellsActive() ) {
return;
}
Opm::DeferredLogger local_deferredLogger;
updatePerforationIntensiveQuantities();
int exception_thrown = 0;
try {
if (iterationIdx == 0) {
calculateExplicitQuantities(local_deferredLogger);
prepareTimeStep(local_deferredLogger);
}
updateWellControls(local_deferredLogger, true);
// only check group controls for iterationIdx smaller then nupcol
const int reportStepIdx = ebosSimulator_.episodeIndex();
const int nupcol = schedule().getNupcol(reportStepIdx);
if (iterationIdx < nupcol) {
if( localWellsActive() ) {
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
std::vector<double> groupTargetReduction(numPhases(), 0.0);
wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, well_state_, groupTargetReduction);
std::vector<double> groupTargetReductionInj(numPhases(), 0.0);
wellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, well_state_, groupTargetReductionInj);
std::vector<double> rein(numPhases(), 0.0);
wellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, well_state_, rein);
double resv = 0.0;
wellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_, resv);
}
}
// Set the well primary variables based on the value of well solutions
initPrimaryVariablesEvaluation();
std::vector< Scalar > B_avg(numComponents(), Scalar() );
computeAverageFormationFactor(B_avg);
if (param_.solve_welleq_initially_ && iterationIdx == 0) {
// solve the well equations as a pre-processing step
last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
if (initial_step_) {
// update the explicit quantities to get the initial fluid distribution in the well correct.
calculateExplicitQuantities(local_deferredLogger);
prepareTimeStep(local_deferredLogger);
last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
initial_step_ = false;
}
// TODO: should we update the explicit related here again, or even prepareTimeStep().
// basically, this is a more updated state from the solveWellEq based on fixed
// reservoir state, will tihs be a better place to inialize the explict information?
}
assembleWellEq(B_avg, dt, local_deferredLogger);
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(local_deferredLogger, exception_thrown, "assemble() failed.", terminal_output_);
last_report_.converged = true;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger)
{
for (auto& well : well_container_) {
well->assembleWellEq(ebosSimulator_, B_avg, dt, well_state_, 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);
}
}
// 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;
int exception_thrown = 0;
try {
if (localWellsActive()) {
for (auto& well : well_container_) {
well->recoverWellSolutionAndUpdateWellState(x, well_state_, local_deferredLogger);
}
}
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(local_deferredLogger, exception_thrown, "recoverWellSolutionAndUpdateWellState() failed.", 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>
SimulatorReport
BlackoilWellModel<TypeTag>::
solveWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger)
{
WellState well_state0 = well_state_;
const int max_iter = param_.max_welleq_iter_;
int it = 0;
bool converged;
int exception_thrown = 0;
do {
try {
assembleWellEq(B_avg, dt, deferred_logger);
} catch (std::exception& e) {
exception_thrown = 1;
}
// We need to check on all processes, as getWellConvergence() below communicates on all processes.
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "solveWellEq() failed.", terminal_output_);
const auto report = getWellConvergence(B_avg);
converged = report.converged();
if (converged) {
break;
}
try {
if( localWellsActive() )
{
for (auto& well : well_container_) {
well->solveEqAndUpdateWellState(well_state_, deferred_logger);
}
}
// updateWellControls uses communication
// Therefore the following is executed if there
// are active wells anywhere in the global domain.
if( wellsActive() )
{
updateWellControls(deferred_logger, /*don't switch group controls*/false);
initPrimaryVariablesEvaluation();
}
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "solveWellEq() failed.", terminal_output_);
++it;
} while (it < max_iter);
try {
if (converged) {
if (terminal_output_) {
deferred_logger.debug("Well equation solution gets converged with " + std::to_string(it) + " iterations");
}
} else {
if (terminal_output_) {
deferred_logger.debug("Well equation solution failed in getting converged with " + std::to_string(it) + " iterations");
}
well_state_ = well_state0;
updatePrimaryVariables(deferred_logger);
}
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "solveWellEq() failed.", terminal_output_);
SimulatorReport report;
report.converged = converged;
report.total_well_iterations = it;
return report;
}
template<typename TypeTag>
ConvergenceReport
BlackoilWellModel<TypeTag>::
getWellConvergence(const std::vector<Scalar>& B_avg) 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(well_state_, 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());
}
}
}
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_, well_state_, deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateWellControls(Opm::DeferredLogger& deferred_logger, const bool checkGroupControl)
{
// 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 ;
// update group controls
if (checkGroupControl) {
const int reportStepIdx = ebosSimulator_.episodeIndex();
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
checkGroupConstraints(fieldGroup, deferred_logger);
}
for (const auto& well : well_container_) {
well->updateWellControl(ebosSimulator_, well_state_, deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
{
Opm::DeferredLogger local_deferredLogger;
for (const auto& well : well_container_) {
well->updateWellTestState(well_state_, simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
}
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_copy = well_state_;
// average B factors are required for the convergence checking of well equations
// Note: this must be done on all processes, even those with
// no wells needing testing, otherwise we will have locking.
std::vector< Scalar > B_avg(numComponents(), Scalar() );
computeAverageFormationFactor(B_avg);
const Opm::SummaryConfig& summaryConfig = ebosSimulator_.vanguard().summaryConfig();
const bool write_restart_file = ebosSimulator_.vanguard().eclState().getRestartConfig().getWriteRestartFile(reportStepIdx);
int exception_thrown = 0;
try {
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());
const Well& eclWell = well->wellEcl();
bool needPotentialsForGuideRate = eclWell.getGuideRatePhase() == Well::GuideRateTarget::UNDEFINED;
if (write_restart_file || needed_for_summary || needPotentialsForGuideRate)
{
std::vector<double> potentials;
well->computeWellPotentials(ebosSimulator_, B_avg, well_state_copy, 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]);
}
}
} // end of for (int w = 0; w < nw; ++w)
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "computeWellPotentials() failed.", terminal_output_);
// Store it in the well state
well_state_.wellPotentials() = well_potentials;
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
prepareTimeStep(Opm::DeferredLogger& deferred_logger)
{
int exception_thrown = 0;
try {
for (const auto& well : well_container_) {
well->checkWellOperability(ebosSimulator_, well_state_, 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 (well_state_.effectiveEventsOccurred(w) ) {
well->updateWellStateWithTarget(ebosSimulator_, well_state_, 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 (well_state_.effectiveEventsOccurred(w) ) {
well_state_.setEffectiveEventsOccurred(w, false);
}
} // end of for (const auto& well : well_container_)
updatePrimaryVariables(deferred_logger);
} catch (std::exception& e) {
exception_thrown = 1;
}
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "prepareTimestep() failed.", terminal_output_);
}
template<typename TypeTag>
const typename BlackoilWellModel<TypeTag>::WellState&
BlackoilWellModel<TypeTag>::
wellState() const { return well_state_; }
template<typename TypeTag>
const typename BlackoilWellModel<TypeTag>::WellState&
BlackoilWellModel<TypeTag>::
wellState(const WellState& well_state OPM_UNUSED) const { return wellState(); }
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) {
for (unsigned i = 0; i < number_of_cells_; ++i) {
cartesian_to_compressed_[global_cell[i]] = i;
}
}
else {
for (unsigned i = 0; i < number_of_cells_; ++i) {
cartesian_to_compressed_[i] = i;
}
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
computeRepRadiusPerfLength(const Grid& grid, Opm::DeferredLogger& deferred_logger)
{
for (const auto& well : well_container_) {
well->computeRepRadiusPerfLength(grid, cartesian_to_compressed_, deferred_logger);
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
computeAverageFormationFactor(std::vector<Scalar>& B_avg) const
{
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 (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_nc_;
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updatePrimaryVariables(Opm::DeferredLogger& deferred_logger)
{
for (const auto& well : well_container_) {
well->updatePrimaryVariables(well_state_, 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 (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] =
grid.cellCenterDepth(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);
}
}
// convert well data from opm-common to well state from opm-core
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
wellsToState( const data::Wells& wells,
const PhaseUsage& phases,
const bool handle_ms_well,
WellStateFullyImplicitBlackoil& state) const
{
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 : state.wellMap() ) {
const auto well_index = wm.second[ 0 ];
const auto& well = wells.at( wm.first );
state.bhp()[ well_index ] = well.bhp;
state.temperature()[ well_index ] = well.temperature;
//state.currentInjectionControls()[ well_index ] = static_cast<Opm::Well::InjectorCMode>(well.injectionControl);
//state.currentProductionControls()[ well_index ] = static_cast<Well::ProducerCMode>(well.productionControl);
const auto wellrate_index = well_index * np;
for( size_t i = 0; i < phs.size(); ++i ) {
assert( well.rates.has( phs[ i ] ) );
state.wellRates()[ wellrate_index + i ] = well.rates.get( phs[ i ] );
}
const auto perforation_pressure = []( const data::Connection& comp ) {
return comp.pressure;
};
const auto perforation_reservoir_rate = []( const data::Connection& comp ) {
return comp.reservoir_rate;
};
std::transform( well.connections.begin(),
well.connections.end(),
state.perfPress().begin() + wm.second[ 1 ],
perforation_pressure );
std::transform( well.connections.begin(),
well.connections.end(),
state.perfRates().begin() + wm.second[ 1 ],
perforation_reservoir_rate );
int local_comp_index = 0;
for (const data::Connection& comp : well.connections) {
const int global_comp_index = wm.second[1] + local_comp_index;
for (int phase_index = 0; phase_index < np; ++phase_index) {
state.perfPhaseRates()[global_comp_index*np + phase_index] = comp.rates.get(phs[phase_index]);
}
++local_comp_index;
}
if (handle_ms_well && !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 = state.topSegmentIndex(well_index);
const auto& segments = 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
state.segPress()[top_segment_index + segment_index] = segment.second.pressure;
const auto& segment_rates = segment.second.rates;
for (int p = 0; p < np; ++p) {
state.segRates()[(top_segment_index + segment_index) * np + p] = segment_rates.get(phs[p]);
}
}
}
}
}
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>::
checkGroupConstraints(const Group& group, Opm::DeferredLogger& deferred_logger) {
// call recursively
const int reportStepIdx = ebosSimulator_.episodeIndex();
for (const std::string& groupName : group.groups()) {
checkGroupConstraints( schedule().getGroup(groupName, reportStepIdx), deferred_logger);
}
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
auto& well_state = well_state_;
if (group.isInjectionGroup())
{
const auto controls = group.injectionControls(summaryState);
int phasePos;
switch (controls.phase) {
case Phase::WATER:
{
phasePos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
break;
}
case Phase::OIL:
{
phasePos = phase_usage_.phase_pos[BlackoilPhases::Liquid];
break;
}
case Phase::GAS:
{
phasePos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
break;
}
default:
throw("Expected WATER, OIL or GAS as type for group injector: " + group.name());
}
if (group.has_control(Group::InjectionCMode::NONE))
{
// do nothing??
}
if (group.has_control(Group::InjectionCMode::RATE))
{
double current_rate = 0.0;
current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
if (controls.surface_max_rate < current_rate) {
actionOnBrokenConstraints(group, Group::InjectionCMode::RATE, reportStepIdx, deferred_logger);
}
}
if (group.has_control(Group::InjectionCMode::RESV))
{
double current_rate = 0.0;
current_rate += wellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
if (controls.resv_max_rate < current_rate) {
actionOnBrokenConstraints(group, Group::InjectionCMode::RESV, reportStepIdx, deferred_logger);
} }
if (group.has_control(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);
double current_rate = 0.0;
current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
if (controls.target_reinj_fraction*production_Rate < current_rate) {
actionOnBrokenConstraints(group, Group::InjectionCMode::REIN, reportStepIdx, deferred_logger);
} }
if (group.has_control(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);
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);
if (controls.target_void_fraction*voidage_rate < total_rate) {
actionOnBrokenConstraints(group, Group::InjectionCMode::VREP, reportStepIdx, deferred_logger);
}
}
if (group.has_control(Group::InjectionCMode::FLD))
{
// do nothing???
//OPM_THROW(std::runtime_error, "Group " + group.name() + "FLD control for injecting groups not implemented" );
}
} else if (group.isProductionGroup())
{
const auto controls = group.productionControls(summaryState);
if (group.has_control(Group::ProductionCMode::NONE))
{
}
if (group.has_control(Group::ProductionCMode::ORAT))
{
double current_rate = 0.0;
current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
if (controls.oil_target < current_rate ) {
actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::ORAT, reportStepIdx, deferred_logger);
}
}
if (group.has_control(Group::ProductionCMode::WRAT))
{
double current_rate = 0.0;
current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
if (controls.water_target < current_rate ) {
actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::WRAT, reportStepIdx, deferred_logger);
}
}
if (group.has_control(Group::ProductionCMode::GRAT))
{
double current_rate = 0.0;
current_rate += wellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
if (controls.gas_target < current_rate ) {
actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::GRAT, reportStepIdx, deferred_logger);
}
}
if (group.has_control(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);
if (controls.liquid_target < current_rate ) {
actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::LRAT, reportStepIdx, deferred_logger);
}
}
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))
{
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);
if (controls.resv_target < current_rate ) {
actionOnBrokenConstraints(group, controls.exceed_action, Group::ProductionCMode::RESV, reportStepIdx, deferred_logger);
}
}
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);
}
if (group.has_control(Group::ProductionCMode::FLD))
{
// do nothing???
}
} else {
//neither production or injecting group FIELD?
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
actionOnBrokenConstraints(const Group& group, const Group::ExceedAction& exceed_action, const Group::ProductionCMode& newControl, const int reportStepIdx, Opm::DeferredLogger& deferred_logger) {
auto& well_state = well_state_;
const Group::ProductionCMode& oldControl = well_state.currentProductionGroupControl(group.name());
std::ostringstream ss;
if (oldControl != newControl) {
const std::string from = Group::ProductionCMode2String(oldControl);
ss << "Group " << group.name() << " exceeding "
<< from << " limit \n";
}
switch(exceed_action) {
case Group::ExceedAction::NONE: {
OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit NONE not implemented", deferred_logger);
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) {
well_state.setCurrentProductionGroupControl(group.name(), newControl);
ss << "Switching control mode for group to " << Group::ProductionCMode2String(newControl);
}
wellGroupHelpers::setGroupControl(group, schedule(), reportStepIdx, false, well_state, ss);
break;
}
default:
throw("Invalid procedure for maximum rate limit selected for group" + group.name());
}
auto cc = Dune::MPIHelper::getCollectiveCommunication();
if (cc.size() > 1) {
ss << " on rank " << cc.rank();
}
if (!ss.str().empty())
deferred_logger.info(ss.str());
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
actionOnBrokenConstraints(const Group& group, const Group::InjectionCMode& newControl, const int reportStepIdx, Opm::DeferredLogger& deferred_logger) {
auto& well_state = well_state_;
const Group::InjectionCMode& oldControl = well_state.currentInjectionGroupControl(group.name());
std::ostringstream ss;
if (oldControl != newControl) {
const std::string from = Group::InjectionCMode2String(oldControl);
ss << "Group " << group.name() << " exceeding "
<< from << " limit \n";
ss << "Switching control mode for group to " << Group::InjectionCMode2String(newControl);
auto cc = Dune::MPIHelper::getCollectiveCommunication();
if (cc.size() > 1) {
ss << " on rank " << cc.rank();
}
well_state.setCurrentInjectionGroupControl(group.name(), newControl);
}
wellGroupHelpers::setGroupControl(group, schedule(), reportStepIdx, /*isInjector*/true, well_state, ss);
if (!ss.str().empty())
deferred_logger.info(ss.str());
}
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;
const Group::InjectionCMode& currentGroupControl = wellState.currentInjectionGroupControl(group.name());
if( currentGroupControl == Group::InjectionCMode::REIN ) {
int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
double gasProductionRate = wellGroupHelpers::sumWellRates(group, schedule, wellState, reportStepIdx, gasPos, /*isInjector*/false);
double solventProductionRate = wellGroupHelpers::sumSolventRates(group, schedule, wellState, reportStepIdx, /*isInjector*/false);
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);
}
}
} // namespace Opm