OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Avoid using the Wells struct.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2019-10-23 09:09:45 +02:00
parent cab0724a26
commit 87188f5862
15 changed files with 782 additions and 702 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -187,6 +187,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/utils/gatherDeferredLogger.hpp
opm/simulators/utils/moduleVersion.hpp
opm/simulators/utils/ParallelRestart.hpp
opm/simulators/wells/PerforationData.hpp
opm/simulators/wells/RateConverter.hpp
opm/simulators/wells/SimFIBODetails.hpp
opm/simulators/wells/WellConnectionAuxiliaryModule.hpp

13
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -40,6 +40,7 @@
#include <opm/core/wells.h>
#include <opm/core/wells/WellCollection.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/wells/VFPInjProperties.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp>
#include <opm/simulators/flow/countGlobalCells.hpp>
@ -233,8 +234,10 @@ namespace Opm {
protected:
Simulator& ebosSimulator_;
std::unique_ptr<WellsManager> wells_manager_;
std::vector< Well > wells_ecl_;
std::vector< std::vector<PerforationData> > well_perf_data_;
std::vector<int> first_perf_index_;
bool wells_active_;
@ -246,8 +249,10 @@ namespace Opm {
std::vector<bool> is_cell_perforated_;
void initializeWellPerfData();
// create the well container
std::vector<WellInterfacePtr > createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger);
std::vector<WellInterfacePtr > createWellContainer(const int time_step);
WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, Opm::DeferredLogger& deferred_logger) const;
@ -278,8 +283,6 @@ namespace Opm {
// used to better efficiency of calcuation
mutable BVector scaleAddRes_;
const Wells* wells() const { return wells_manager_->c_wells(); }
const Grid& grid() const
{ return ebosSimulator_.vanguard().grid(); }
@ -373,7 +376,7 @@ namespace Opm {
WellStateFullyImplicitBlackoil& state ) const;
// whether there exists any multisegment well open on this process
bool anyMSWellOpenLocal(const Wells* wells) const;
bool anyMSWellOpenLocal() const;
const Well& getWellEcl(const std::string& well_name) const;

320
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -37,6 +37,19 @@ namespace Opm {
// 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>
@ -46,31 +59,15 @@ namespace Opm {
{
const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState();
gravity_ = ebosSimulator_.problem().gravity()[2];
extractLegacyCellPvtRegionIndex_();
extractLegacyDepth_();
phase_usage_ = phaseUsageFromDeck(eclState);
const auto& gridView = ebosSimulator_.gridView();
// 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()
number_of_cells_ = gridView.size(/*codim=*/0);
global_nc_ = gridView.comm().sum(number_of_cells_);
gravity_ = ebosSimulator_.problem().gravity()[2];
extractLegacyCellPvtRegionIndex_();
extractLegacyDepth_();
initial_step_ = true;
const auto& grid = ebosSimulator_.vanguard().grid();
const auto& cartDims = Opm::UgGridHelpers::cartDims(grid);
setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid),
cartDims[0]*cartDims[1]*cartDims[2]);
// add the eWoms auxiliary module for the wells to the list
ebosSimulator_.model().addAuxiliaryModule(this);
@ -210,24 +207,18 @@ namespace Opm {
Opm::DeferredLogger local_deferredLogger;
const Grid& grid = ebosSimulator_.vanguard().grid();
const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
const auto& eclState = ebosSimulator_.vanguard().eclState();
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
wells_ecl_ = schedule().getWells(timeStepIdx);
// Create wells and well state.
wells_manager_.reset( new WellsManager (eclState,
schedule(),
summaryState,
timeStepIdx,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::dimensions(grid),
Opm::UgGridHelpers::cell2Faces(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
grid.comm().size() > 1,
defunct_well_names) );
// 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.
@ -269,35 +260,22 @@ namespace Opm {
}
cellPressures[cellIdx] = perf_pressure;
}
well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_);
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(wells())) { // if we use MultisegmentWell model
well_state_.initWellStateMSWell(wells(), wells_ecl_, phase_usage_, &previous_well_state_);
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()->number_of_wells;
const int nw = wells_ecl_.size();
for (int w = 0; w <nw; ++w) {
const int nw_wells_ecl = wells_ecl_.size();
int index_well_ecl = 0;
const std::string well_name(wells()->name[w]);
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_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
}
const auto well = wells_ecl_[index_well_ecl];
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))
if(!schedule().hasWellEvent(well.name(), effective_events_mask, timeStepIdx))
continue;
if (well.isProducer()) {
@ -347,7 +325,7 @@ namespace Opm {
wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
// create the well container
well_container_ = createWellContainer(reportStepIdx, local_deferredLogger);
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
@ -400,9 +378,16 @@ namespace Opm {
//compute well guideRates
const int np = numPhases();
for (const auto& well : well_container_) {
const double& oilpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Liquid]];
const double& gaspot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Vapour]];
const double& waterpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Aqua]];
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);
@ -470,7 +455,7 @@ namespace Opm {
Opm::DeferredLogger local_deferredLogger;
for (const auto& well : well_container_) {
if (GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well->wellType() == INJECTOR) {
if (GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well->isInjector()) {
well->updateWaterThroughput(dt, well_state_);
}
}
@ -538,8 +523,6 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
initFromRestartFile(const RestartValue& restartValues)
{
const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
// 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
@ -547,29 +530,25 @@ namespace Opm {
const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0);
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
WellsManager wellsmanager(eclState(),
schedule(),
summaryState,
report_step,
Opm::UgGridHelpers::numCells(grid()),
Opm::UgGridHelpers::globalCell(grid()),
Opm::UgGridHelpers::cartDims(grid()),
Opm::UgGridHelpers::dimensions(grid()),
Opm::UgGridHelpers::cell2Faces(grid()),
Opm::UgGridHelpers::beginFaceCentroids(grid()),
grid().comm().size() > 1,
defunctWellNames);
// 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);
}
const Wells* wells = wellsmanager.c_wells();
initializeWellPerfData();
wells_ecl_ = schedule().getWells(report_step);
const int nw = wells->number_of_wells;
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(wells));
well_state_.resize(wells, wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage); // Resize for restart step
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_);
}
@ -578,20 +557,7 @@ namespace Opm {
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 int nw_wells_ecl = wells_ecl_.size();
int index_well_ecl = 0;
const std::string well_name(wells->name[w]);
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_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
}
const auto& well = wells_ecl_[index_well_ecl];
const auto& well = wells_ecl_[w];
if (well.isProducer()) {
const auto controls = well.productionControls(summaryState);
@ -613,10 +579,59 @@ namespace Opm {
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, Opm::DeferredLogger& deferred_logger)
createWellContainer(const int time_step)
{
std::vector<WellInterfacePtr> well_container;
@ -624,27 +639,9 @@ namespace Opm {
if (nw > 0) {
well_container.reserve(nw);
// With the following way, it will have the same order with wells struct
// Hopefully, it can generate the same residual history with master branch
for (int w = 0; w < nw; ++w) {
const std::string well_name = std::string(wells()->name[w]);
// finding the location of the well in wells_ecl
const int nw_wells_ecl = wells_ecl_.size();
int index_well = 0;
for (; index_well < nw_wells_ecl; ++index_well) {
if (well_name == wells_ecl_[index_well].name()) {
break;
}
}
// It should be able to find in wells_ecl.
if (index_well == nw_wells_ecl) {
OPM_DEFLOG_THROW(std::runtime_error, "Could not find well " + well_name + " in wells_ecl ", deferred_logger);
}
const Well& well_ecl = wells_ecl_[index_well];
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)) {
@ -663,7 +660,6 @@ namespace Opm {
}
}
// TODO: should we do this for all kinds of closing reasons?
// something like wellTestState_.hasWell(well_name)?
bool wellIsStopped = false;
@ -688,15 +684,31 @@ namespace Opm {
}
// Use the pvtRegionIdx from the top cell
const int well_cell_top = wells()->well_cells[wells()->well_connpos[w]];
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, wells(),
param_, *rateConverter_, pvtreg, numComponents() ) );
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, wells(),
param_, *rateConverter_, pvtreg, numComponents() ) );
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();
@ -732,30 +744,32 @@ namespace Opm {
const Well& well_ecl = wells_ecl_[index_well_ecl];
// Finding the location of the well in wells struct.
const int nw = numLocalWells();
int well_index_wells = -999;
for (int w = 0; w < nw; ++w) {
if (well_name == std::string(wells()->name[w])) {
well_index_wells = w;
break;
}
}
if (well_index_wells < 0) {
OPM_DEFLOG_THROW(std::logic_error, "Could not find the well " << well_name << " in the well struct ", deferred_logger);
}
// Use the pvtRegionIdx from the top cell
const int well_cell_top = wells()->well_cells[wells()->well_connpos[well_index_wells]];
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, wells(),
param_, *rateConverter_, pvtreg, numComponents() ) );
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, wells(),
param_, *rateConverter_, pvtreg, numComponents() ) );
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]));
}
}
@ -1183,10 +1197,10 @@ namespace Opm {
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->wellType() == INJECTOR) ||
summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->isInjector()) ||
((summaryConfig.hasSummaryKey( "WWPP:" + well->name()) ||
summaryConfig.hasSummaryKey( "WOPP:" + well->name()) ||
summaryConfig.hasSummaryKey( "WGPP:" + well->name())) && well->wellType() == PRODUCER);
summaryConfig.hasSummaryKey( "WGPP:" + well->name())) && well->isProducer());
const Well& eclWell = well->wellEcl();
bool needPotentialsForGuideRate = eclWell.getGuideRatePhase() == Well::GuideRateTarget::UNDEFINED;
@ -1411,14 +1425,14 @@ namespace Opm {
int
BlackoilWellModel<TypeTag>:: numLocalWells() const
{
return wells() ? wells()->number_of_wells : 0;
return wells_ecl_.size();
}
template<typename TypeTag>
int
BlackoilWellModel<TypeTag>:: numPhases() const
BlackoilWellModel<TypeTag>::numPhases() const
{
return wells() ? wells()->number_of_phases : 1;
return phase_usage_.num_phases;
}
template<typename TypeTag>
@ -1560,22 +1574,14 @@ namespace Opm {
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
anyMSWellOpenLocal(const Wells* wells) const
anyMSWellOpenLocal() const
{
bool any_ms_well_open = false;
const int nw = wells->number_of_wells;
for (int w = 0; w < nw; ++w) {
const std::string well_name = std::string(wells->name[w]);
const Well& well_ecl = getWellEcl(well_name);
if (well_ecl.isMultiSegment() ) {
any_ms_well_open = true;
break;
for (const auto& well : wells_ecl_) {
if (well.isMultiSegment()) {
return true;
}
}
return any_ms_well_open;
return false;
}

10
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -98,11 +98,15 @@ namespace Opm
// TODO: for now, we only use one type to save some implementation efforts, while improve later.
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
MultisegmentWell(const Well& well, const int time_step, const Wells* wells,
MultisegmentWell(const Well& well, const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components);
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
@ -185,7 +189,6 @@ namespace Opm
using Base::well_cells_;
using Base::param_;
using Base::well_index_;
using Base::well_type_;
using Base::first_perf_;
using Base::saturation_table_number_;
using Base::well_efficiency_factor_;
@ -355,7 +358,6 @@ namespace Opm
void assembleGroupProductionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
void assembleGroupInjectionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, const Well::InjectorType& injectorType, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
void assemblePressureEq(const int seg) const;
// hytrostatic pressure loss

60
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -28,12 +28,16 @@ namespace Opm
template <typename TypeTag>
MultisegmentWell<TypeTag>::
MultisegmentWell(const Well& well, const int time_step, const Wells* wells,
MultisegmentWell(const Well& well, const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components)
: Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components)
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data)
: Base(well, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, first_perf_index, perf_data)
, segment_perforations_(numberOfSegments())
, segment_inlets_(numberOfSegments())
, cell_perforation_depth_diffs_(number_of_perforations_, 0.0)
@ -756,7 +760,7 @@ namespace Opm
primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
}
} else { // total_seg_rate == 0
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
// only single phase injection handled
auto phase = well.getInjectionProperties().injectorType;
@ -776,7 +780,7 @@ namespace Opm
}
}
} else if (well_type_ == PRODUCER) { // producers
} else if (this->isProducer()) { // producers
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
}
@ -942,7 +946,7 @@ namespace Opm
// make sure that no injector produce and no producer inject
if (seg == 0) {
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
primary_variables_[seg][GTotal] = std::max( primary_variables_[seg][GTotal], 0.0);
} else {
primary_variables_[seg][GTotal] = std::min( primary_variables_[seg][GTotal], 0.0);
@ -1186,7 +1190,7 @@ namespace Opm
// producing perforations
if ( drawdown > 0.0) {
// Do nothing is crossflow is not allowed
if (!allow_cf && well_type_ == INJECTOR) {
if (!allow_cf && this->isInjector()) {
return;
}
@ -1206,7 +1210,7 @@ namespace Opm
}
} else { // injecting perforations
// Do nothing if crossflow is not allowed
if (!allow_cf && well_type_ == PRODUCER) {
if (!allow_cf && this->isProducer()) {
return;
}
@ -1263,7 +1267,7 @@ namespace Opm
} // end for injection perforations
// calculating the perforation solution gas rate and solution oil rates
if (well_type_ == PRODUCER) {
if (this->isProducer()) {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
@ -1516,7 +1520,7 @@ namespace Opm
// the result will contain the derivative with resepct to GTotal in segment seg,
// and the derivatives with respect to WFrac GFrac in segment seg_upwind.
// the derivative with respect to SPres should be zero.
if (seg == 0 && well_type_ == INJECTOR) {
if (seg == 0 && this->isInjector()) {
const Well& well = Base::wellEcl();
auto phase = well.getInjectionProperties().injectorType;
@ -1630,7 +1634,7 @@ namespace Opm
if (wellIsStopped_) {
control_eq = getSegmentGTotal(0);
} else if (well.isInjector() ) {
} else if (this->isInjector() ) {
const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
const auto controls = well.injectionControls(summaryState);
@ -1711,13 +1715,16 @@ namespace Opm
rates[ Gas ] = getSegmentRate(0, Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
}
control_eq = getSegmentPressure(0) - calculateBhpFromThp(rates, well, summaryState, deferred_logger);
break; }
break;
}
case Well::InjectorCMode::BHP:
{
const auto& bhp = controls.bhp_limit;
control_eq = getSegmentPressure(0) - bhp;
break;
}
case Well::InjectorCMode::GRUP:
{
assert(well.isAvailableForGroupControl());
@ -1725,6 +1732,7 @@ namespace Opm
assembleGroupInjectionControl(group, well_state, schedule, summaryState, controls.injector_type, control_eq, efficiencyFactor, deferred_logger);
break;
}
case Well::InjectorCMode::CMODE_UNDEFINED:
{
OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + name(), deferred_logger);
@ -1823,7 +1831,8 @@ namespace Opm
rates[ Gas ] = getSegmentRate(0, Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx));
}
control_eq = getSegmentPressure(0) - calculateBhpFromThp(rates, well, summaryState, deferred_logger);
break; }
break;
}
case Well::ProducerCMode::GRUP:
{
assert(well.isAvailableForGroupControl());
@ -1844,7 +1853,6 @@ namespace Opm
}
}
// using control_eq to update the matrix and residuals
resWell_[0][SPres] = control_eq.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
@ -1903,14 +1911,14 @@ namespace Opm
const double rho = segment_densities_[0].value();
double thp = 0.0;
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
const int table_id = well_ecl_.vfp_table_number();
const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (well_type_ == PRODUCER) {
else if (this->isProducer()) {
const int table_id = well_ecl_.vfp_table_number();
const double alq = well_ecl_.alq_value();
const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
@ -2675,7 +2683,7 @@ namespace Opm
computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
// updating the solution gas rate and solution oil rate
if (well_type_ == PRODUCER) {
if (this->isProducer()) {
well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate;
well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate;
}
@ -2762,8 +2770,8 @@ namespace Opm
// for now, if there is one perforation can produce/inject in the correct
// direction, we consider this well can still produce/inject.
// TODO: it can be more complicated than this to cause wrong-signed rates
if ( (drawdown < 0. && well_type_ == INJECTOR) ||
(drawdown > 0. && well_type_ == PRODUCER) ) {
if ( (drawdown < 0. && this->isInjector()) ||
(drawdown > 0. && this->isProducer()) ) {
all_drawdown_wrong_direction = false;
break;
}
@ -3056,9 +3064,8 @@ namespace Opm
{
double control_tolerance = 0.;
const auto& well = well_ecl_;
const int well_index = index_of_well_;
if (well.isInjector() )
if (this->isInjector() )
{
const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
switch(current) {
@ -3080,7 +3087,7 @@ namespace Opm
}
}
if (well.isProducer() )
if (this->isProducer() )
{
const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
switch(current) {
@ -3124,9 +3131,8 @@ namespace Opm
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
const auto& well = well_ecl_;
const int well_index = index_of_well_;
if (well.isInjector() )
if (this->isInjector() )
{
const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
switch(current) {
@ -3152,7 +3158,7 @@ namespace Opm
}
}
if (well.isProducer() )
if (this->isProducer() )
{
const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
switch(current) {
@ -3208,8 +3214,8 @@ namespace Opm
// special treatment is needed for segment 0
if (seg == 0) {
// we are not supposed to have injecting producers and producing injectors
assert( ! (well_type_ == PRODUCER && primary_variables_evaluation_[seg][GTotal] > 0.) );
assert( ! (well_type_ == INJECTOR && primary_variables_evaluation_[seg][GTotal] < 0.) );
assert( ! (this->isProducer() && primary_variables_evaluation_[seg][GTotal] > 0.) );
assert( ! (this->isInjector() && primary_variables_evaluation_[seg][GTotal] < 0.) );
upwinding_segments_[seg] = seg;
continue;
}

36
opm/simulators/wells/PerforationData.hpp Normal file
View File

@ -0,0 +1,36 @@
/*
Copyright 2019 SINTEF Digital, Mathematics and Cybernetics.
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/>.
*/
#ifndef OPM_PERFORATIONDATA_HEADER_INCLUDED
#define OPM_PERFORATIONDATA_HEADER_INCLUDED
namespace Opm
{
/// Static data associated with a well perforation.
struct PerforationData
{
int cell_index;
double connection_transmissibility_factor;
int satnum_id;
};
} // namespace Opm
#endif // OPM_PERFORATIONDATA_HEADER_INCLUDED

11
opm/simulators/wells/StandardWell.hpp
View File

@ -140,11 +140,15 @@ namespace Opm
using Base::contiFoamEqIdx;
static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
StandardWell(const Well& well, const int time_step, const Wells* wells,
StandardWell(const Well& well, const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components);
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
@ -231,10 +235,8 @@ namespace Opm
using Base::number_of_perforations_;
using Base::number_of_phases_;
using Base::saturation_table_number_;
using Base::comp_frac_;
using Base::well_index_;
using Base::index_of_well_;
using Base::well_type_;
using Base::num_components_;
using Base::connectionRates_;
@ -398,7 +400,6 @@ namespace Opm
void assembleGroupProductionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
void assembleGroupInjectionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, const Well::InjectorType& injectorType, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
// handle the non reasonable fractions due to numerical overshoot
void processFractions() const;

216
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -28,12 +28,16 @@ namespace Opm
template<typename TypeTag>
StandardWell<TypeTag>::
StandardWell(const Well& well, const int time_step, const Wells* wells,
StandardWell(const Well& well, const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components)
: Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components)
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data)
: Base(well, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, first_perf_index, perf_data)
, perf_densities_(number_of_perforations_)
, perf_pressure_diffs_(number_of_perforations_)
, F0_(numWellConservationEq)
@ -68,7 +72,7 @@ namespace Opm
}
// counting/updating primary variable numbers
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
// adding a primary variable for water perforation rate per connection
numWellEq_ += number_of_perforations_;
// adding a primary variable for skin pressure per connection
@ -186,25 +190,33 @@ namespace Opm
// Note: currently, the WQTotal definition is still depends on Injector/Producer.
assert(comp_idx < num_components_);
if (well_type_ == INJECTOR) { // only single phase injection
// TODO: using comp_frac here is dangerous, it should be changed later
// Most likely, it should be changed to use distr, or at least, we need to update comp_frac_ based on distr
// while solvent might complicate the situation
const auto pu = phaseUsage();
const int legacyCompIdx = ebosCompIdxToFlowCompIdx(comp_idx);
double comp_frac = 0.0;
if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent
comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
} else if (legacyCompIdx == pu.phase_pos[ Gas ]) {
comp_frac = comp_frac_[legacyCompIdx];
if (has_solvent) {
comp_frac *= (1.0 - wsolvent());
if (this->isInjector()) { // only single phase injection
double inj_frac = 0.0;
switch (this->wellEcl().injectorType()) {
case Well::InjectorType::WATER:
if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx))) {
inj_frac = 1.0;
}
} else {
comp_frac = comp_frac_[legacyCompIdx];
break;
case Well::InjectorType::GAS:
if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent
inj_frac = wsolvent();
} else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) {
inj_frac = has_solvent ? 1.0 - wsolvent() : 1.0;
}
break;
case Well::InjectorType::OIL:
if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) {
inj_frac = 1.0;
}
break;
case Well::InjectorType::MULTI:
// Not supported.
// deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
// "Multi phase injectors are not supported, requested for well " + name());
break;
}
return comp_frac * primary_variables_evaluation_[WQTotal];
return inj_frac * primary_variables_evaluation_[WQTotal];
} else { // producers
return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx);
}
@ -364,7 +376,7 @@ namespace Opm
const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
EvalWell drawdown = pressure - well_pressure;
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index];
drawdown += skin_pressure;
@ -373,7 +385,7 @@ namespace Opm
// producing perforations
if ( drawdown.value() > 0 ) {
//Do nothing if crossflow is not allowed
if (!allow_cf && well_type_ == INJECTOR) {
if (!allow_cf && this->isInjector()) {
return;
}
@ -395,7 +407,7 @@ namespace Opm
cq_s[oilCompIdx] += vap_oil;
// recording the perforation solution gas rate and solution oil rates
if (well_type_ == PRODUCER) {
if (this->isProducer()) {
perf_dis_gas_rate = dis_gas.value();
perf_vap_oil_rate = vap_oil.value();
}
@ -403,7 +415,7 @@ namespace Opm
} else {
//Do nothing if crossflow is not allowed
if (!allow_cf && well_type_ == PRODUCER) {
if (!allow_cf && this->isProducer()) {
return;
}
@ -471,7 +483,7 @@ namespace Opm
}
// calculating the perforation solution gas rate and solution oil rates
if (well_type_ == PRODUCER) {
if (this->isProducer()) {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
@ -553,12 +565,12 @@ namespace Opm
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
// better way to do here is that use the cq_s and then replace the cq_s_water here?
if (has_polymer && this->has_polymermw && well_type_ == INJECTOR) {
if (has_polymer && this->has_polymermw && this->isInjector()) {
handleInjectivityRateAndEquations(intQuants, well_state, perf, cq_s, deferred_logger);
}
// updating the solution gas rate and solution oil rate
if (well_type_ == PRODUCER) {
if (this->isProducer()) {
well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate;
well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate;
}
@ -629,7 +641,7 @@ namespace Opm
}
// change temperature for injecting fluids
if (well_type_ == INJECTOR && cq_s[activeCompIdx] > 0.0){
if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
// only handles single phase injection now
assert(this->well_ecl_.injectorType() != Well::InjectorType::MULTI);
fs.setTemperature(this->well_ecl_.temperature());
@ -655,7 +667,7 @@ namespace Opm
// TODO: the application of well efficiency factor has not been tested with an example yet
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
EvalWell cq_s_poly = cq_s[waterCompIdx] * well_efficiency_factor_;
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
cq_s_poly *= wpolymer();
} else {
cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
@ -671,7 +683,7 @@ namespace Opm
// TODO: the application of well efficiency factor has not been tested with an example yet
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
EvalWell cq_s_foam = cq_s[gasCompIdx] * well_efficiency_factor_;
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
cq_s_foam *= wfoam();
} else {
cq_s_foam *= extendEval(intQuants.foamConcentration());
@ -751,7 +763,7 @@ namespace Opm
if (wellIsStopped_) {
control_eq = getWQTotal();
} else if (well.isInjector()) {
} else if (this->isInjector()) {
const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
const auto& controls = well.injectionControls(summaryState);
switch(current) {
@ -1309,7 +1321,7 @@ namespace Opm
// for injectors, very typical one of the fractions will be one, and it is easy to get zero value
// fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
const double relaxation_factor_fractions = (well_type_ == PRODUCER) ?
const double relaxation_factor_fractions = (this->isProducer()) ?
relaxationFactorFractionsProducer(old_primary_variables, dwells)
: 1.0;
@ -1368,7 +1380,7 @@ namespace Opm
updateExtraPrimaryVariables(const BVectorWell& dwells) const
{
// for the water velocity and skin pressure
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const int wat_vel_index = Bhp + 1 + perf;
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
@ -1520,24 +1532,37 @@ namespace Opm
// calculate the phase rates based on the primary variables
// for producers, this is not a problem, while not sure for injectors here
if (well_type_ == PRODUCER) {
if (this->isProducer()) {
const double g_total = primary_variables_[WQTotal];
for (int p = 0; p < number_of_phases_; ++p) {
well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = g_total * F[p];
}
} else { // injectors
// TODO: using comp_frac_ here is very dangerous, since we do not update it based on the injection phase
// Either we use distr (might conflict with RESV related) or we update comp_frac_ based on the injection phase
for (int p = 0; p < number_of_phases_; ++p) {
const double comp_frac = comp_frac_[p];
well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = comp_frac * primary_variables_[WQTotal];
well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = 0.0;
}
switch (this->wellEcl().injectorType()) {
case Well::InjectorType::WATER:
well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[Water]] = primary_variables_[WQTotal];
break;
case Well::InjectorType::GAS:
well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[Gas]] = primary_variables_[WQTotal];
break;
case Well::InjectorType::OIL:
well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[Oil]] = primary_variables_[WQTotal];
break;
case Well::InjectorType::MULTI:
// Not supported.
deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
"Multi phase injectors are not supported, requested for well " + name());
break;
}
}
updateThp(well_state, deferred_logger);
// other primary variables related to polymer injectivity study
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf];
well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf];
@ -1607,7 +1632,7 @@ namespace Opm
return;
}
if (well.isInjector() )
if (this->isInjector() )
{
const auto& controls = well.injectionControls(summaryState);
@ -1824,7 +1849,7 @@ namespace Opm
// TODO: it only handles the producers for now
// the formular for the injectors are not formulated yet
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
return;
}
@ -1926,7 +1951,7 @@ namespace Opm
}
// focusing on PRODUCER for now
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
return;
}
@ -2105,8 +2130,8 @@ namespace Opm
// for now, if there is one perforation can produce/inject in the correct
// direction, we consider this well can still produce/inject.
// TODO: it can be more complicated than this to cause wrong-signed rates
if ( (drawdown < 0. && well_type_ == INJECTOR) ||
(drawdown > 0. && well_type_ == PRODUCER) ) {
if ( (drawdown < 0. && this->isInjector()) ||
(drawdown > 0. && this->isProducer()) ) {
all_drawdown_wrong_direction = false;
break;
}
@ -2129,15 +2154,15 @@ namespace Opm
std::vector<double> well_rates;
computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
const double sign = (well_type_ == PRODUCER) ? -1. : 1.;
const double sign = (this->isProducer()) ? -1. : 1.;
const double threshold = sign * std::numeric_limits<double>::min();
bool can_produce_inject = false;
for (const auto value : well_rates) {
if (well_type_ == PRODUCER && value < threshold) {
if (this->isProducer() && value < threshold) {
can_produce_inject = true;
break;
} else if (well_type_ == INJECTOR && value > threshold) {
} else if (this->isInjector() && value > threshold) {
can_produce_inject = true;
break;
}
@ -2293,7 +2318,6 @@ namespace Opm
const std::vector<double>& surf_dens_perf)
{
// Verify that we have consistent input.
const int np = number_of_phases_;
const int nperf = number_of_perforations_;
const int num_comp = num_components_;
@ -2336,13 +2360,44 @@ namespace Opm
mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate);
}
} else {
std::fill(mix.begin(), mix.end(), 0.0);
// No flow => use well specified fractions for mix.
for (int component = 0; component < num_comp; ++component) {
if (component < np) {
mix[component] = comp_frac_[ ebosCompIdxToFlowCompIdx(component)];
if (this->isInjector()) {
switch (this->wellEcl().injectorType()) {
case Well::InjectorType::WATER:
mix[FluidSystem::waterCompIdx] = 1.0;
break;
case Well::InjectorType::GAS:
mix[FluidSystem::gasCompIdx] = 1.0;
break;
case Well::InjectorType::OIL:
mix[FluidSystem::oilCompIdx] = 1.0;
break;
case Well::InjectorType::MULTI:
// Not supported.
// deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
// "Multi phase injectors are not supported, requested for well " + name());
break;
}
} else {
assert(this->isProducer());
// Using the preferred phase to decide the mix initialization.
switch (this->wellEcl().getPreferredPhase()) {
case Phase::OIL:
mix[FluidSystem::oilCompIdx] = 1.0;
break;
case Phase::GAS:
mix[FluidSystem::gasCompIdx] = 1.0;
break;
case Phase::WATER:
mix[FluidSystem::waterCompIdx] = 1.0;
break;
default:
// No others supported.
break;
}
}
// intialize 0.0 for comIdx >= np;
}
// Compute volume ratio.
x = mix;
@ -2848,6 +2903,7 @@ namespace Opm
const int well_index = index_of_well_;
const int np = number_of_phases_;
const auto& pu = phaseUsage();
// the weighted total well rate
double total_well_rate = 0.0;
@ -2857,12 +2913,22 @@ namespace Opm
// Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate
// under surface condition is used here
if (well_type_ == INJECTOR) {
primary_variables_[WQTotal] = 0.;
for (int p = 0; p < np; ++p) {
// TODO: the use of comp_frac_ here is dangerous, since the injection phase can be different from
// prefered phasse in WELSPECS, while comp_frac_ only reflect the one specified in WELSPECS
primary_variables_[WQTotal] += well_state.wellRates()[np * well_index + p] * comp_frac_[p];
if (this->isInjector()) {
switch (this->wellEcl().injectorType()) {
case Well::InjectorType::WATER:
primary_variables_[WQTotal] = well_state.wellRates()[np * well_index + pu.phase_pos[Water]];
break;
case Well::InjectorType::GAS:
primary_variables_[WQTotal] = well_state.wellRates()[np * well_index + pu.phase_pos[Gas]];
break;
case Well::InjectorType::OIL:
primary_variables_[WQTotal] = well_state.wellRates()[np * well_index + pu.phase_pos[Oil]];
break;
case Well::InjectorType::MULTI:
// Not supported.
deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED",
"Multi phase injectors are not supported, requested for well " + name());
break;
}
} else {
for (int p = 0; p < np; ++p) {
@ -2871,7 +2937,6 @@ namespace Opm
}
if (std::abs(total_well_rate) > 0.) {
const auto pu = phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate;
}
@ -2882,7 +2947,7 @@ namespace Opm
primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ;
}
} else { // total_well_rate == 0
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
auto phase = well_ecl_.getInjectionProperties().injectorType;
// only single phase injection handled
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
@ -2907,7 +2972,7 @@ namespace Opm
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_type_ == PRODUCER) { // producers
} else if (this->isProducer()) { // producers
// TODO: the following are not addressed for the solvent case yet
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
primary_variables_[WFrac] = 1.0 / np;
@ -2925,7 +2990,7 @@ namespace Opm
primary_variables_[Bhp] = well_state.bhp()[index_of_well_];
// other primary variables related to polymer injection
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf];
primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf];
@ -2965,14 +3030,14 @@ namespace Opm
// TODO: it is possible it should be a Evaluation
const double rho = perf_densities_[0];
if (well.isInjector() )
if (this->isInjector() )
{
const auto& controls = well.injectionControls(summaryState);
const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit) - dp;
}
else if (well.isProducer()) {
else if (this->isProducer()) {
const auto& controls = well.productionControls(summaryState);
const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
@ -3005,14 +3070,14 @@ namespace Opm
const double rho = perf_densities_[0];
double thp = 0.0;
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
const int table_id = well_ecl_.vfp_table_number();
const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (well_type_ == PRODUCER) {
else if (this->isProducer()) {
const int table_id = well_ecl_.vfp_table_number();
const double alq = well_ecl_.alq_value();
const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
@ -3052,7 +3117,7 @@ namespace Opm
// TODO: not sure should based on the well type or injecting/producing peforations
// it can be different for crossflow
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
// assume fully mixing within injecting wellbore
const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
@ -3062,7 +3127,7 @@ namespace Opm
if (PolymerModule::hasPlyshlog()) {
// we do not calculate the shear effects for injection wells when they do not
// inject polymer.
if (well_type_ == INJECTOR && wpolymer() == 0.) {
if (this->isInjector() && wpolymer() == 0.) {
return;
}
// compute the well water velocity with out shear effects.
@ -3398,7 +3463,7 @@ namespace Opm
StandardWell<TypeTag>::
updateWaterThroughput(const double dt, WellState &well_state) const
{
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double perf_water_vel = primary_variables_[Bhp + 1 + perf];
// we do not consider the formation damage due to water flowing from reservoir into wellbore
@ -3477,13 +3542,12 @@ namespace Opm
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
const auto& well = well_ecl_;
const int well_index = index_of_well_;
if (wellIsStopped_) {
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; // use smaller tolerance for zero control?
}
else if (well.isInjector() )
else if (this->isInjector() )
{
const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index];
switch(current) {
@ -3508,7 +3572,7 @@ namespace Opm
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
}
}
else if (well.isProducer() )
else if (this->isProducer() )
{
const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index];
switch(current) {
@ -3563,7 +3627,7 @@ namespace Opm
// if different types of extra equations are involved, this function needs to be refactored further
// checking the convergence of the extra equations related to polymer injectivity
if (this->has_polymermw && well_type_ == INJECTOR) {
if (this->has_polymermw && this->isInjector()) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const int dummy_component = -1;
@ -3612,7 +3676,7 @@ namespace Opm
{
// the source term related to transport of molecular weight
EvalWell cq_s_polymw = cq_s_poly;
if (well_type_ == INJECTOR) {
if (this->isInjector()) {
const int wat_vel_index = Bhp + 1 + perf;
const EvalWell water_velocity = primary_variables_evaluation_[wat_vel_index];
if (water_velocity > 0.) { // injecting
@ -3624,7 +3688,7 @@ namespace Opm
// going-back to the wellbore through injector
cq_s_polymw *= 0.;
}
} else if (well_type_ == PRODUCER) {
} else if (this->isProducer()) {
if (cq_s_polymw < 0.) {
cq_s_polymw *= extendEval(int_quants.polymerMoleWeight() );
} else {

47
opm/simulators/wells/WellInterface.hpp
View File

@ -118,11 +118,15 @@ namespace Opm
compositionSwitchEnabled,
Indices::numPhases >;
/// Constructor
WellInterface(const Well& well, const int time_step, const Wells* wells,
WellInterface(const Well& well, const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components);
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
/// Virutal destructor
virtual ~WellInterface() {}
@ -130,15 +134,18 @@ namespace Opm
/// Well name.
const std::string& name() const;
/// True if the well is an injector.
bool isInjector() const;
/// True if the well is a producer.
bool isProducer() const;
/// Index of well in the wells struct and wellState
int indexOfWell() const;
/// Well cells.
const std::vector<int>& cells() const {return well_cells_; }
/// Well type, INJECTOR or PRODUCER.
WellType wellType() const;
/// Well controls
WellControls* wellControls() const;
@ -286,30 +293,12 @@ namespace Opm
const int current_step_;
// the index of well in Wells struct
int index_of_well_;
// simulation parameters
const ModelParameters& param_;
// well type
// INJECTOR or PRODUCER
enum WellType well_type_;
// number of phases
int number_of_phases_;
// component fractions for each well
// typically, it should apply to injection wells
std::vector<double> comp_frac_;
// number of the perforations for this well
int number_of_perforations_;
// record the index of the first perforation
// of states of individual well.
int first_perf_;
// well index for each perforation
std::vector<double> well_index_;
@ -372,6 +361,18 @@ namespace Opm
const int num_components_;
// number of phases
int number_of_phases_;
// the index of well in Wells struct
int index_of_well_;
// record the index of the first perforation
// of states of individual well.
int first_perf_;
std::vector<int> originalConnectionIndex_;
std::vector<RateVector> connectionRates_;
bool wellIsStopped_;

113
opm/simulators/wells/WellInterface_impl.hpp
View File

@ -27,75 +27,58 @@ namespace Opm
template<typename TypeTag>
WellInterface<TypeTag>::
WellInterface(const Well& well, const int time_step, const Wells* wells,
WellInterface(const Well& well,
const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components)
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data)
: well_ecl_(well)
, current_step_(time_step)
, param_(param)
, rateConverter_(rate_converter)
, pvtRegionIdx_(pvtRegionIdx)
, num_components_(num_components)
, number_of_phases_(num_phases)
, index_of_well_(index_of_well)
, first_perf_(first_perf_index)
{
if (time_step < 0) {
OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
}
if (!wells) {
OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface");
}
ref_depth_ = well.getRefDepth();
const std::string& well_name = well.name();
// looking for the location of the well in the wells struct
int index_well;
for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
if (well_name == std::string(wells->name[index_well])) {
break;
}
}
// should not enter the constructor if the well does not exist in the wells struct
// here, just another assertion.
assert(index_well != wells->number_of_wells);
index_of_well_ = index_well;
well_type_ = wells->type[index_well];
number_of_phases_ = wells->number_of_phases;
// copying the comp_frac
{
comp_frac_.resize(number_of_phases_);
const int index_begin = index_well * number_of_phases_;
std::copy(wells->comp_frac + index_begin,
wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() );
}
ref_depth_ = wells->depth_ref[index_well];
// We do not want to count SHUT perforations here, so
// it would be wrong to use wells.getConnections().size().
number_of_perforations_ = perf_data.size();
// perforations related
{
const int perf_index_begin = wells->well_connpos[index_well];
const int perf_index_end = wells->well_connpos[index_well + 1];
number_of_perforations_ = perf_index_end - perf_index_begin;
first_perf_ = perf_index_begin;
well_cells_.resize(number_of_perforations_);
std::copy(wells->well_cells + perf_index_begin,
wells->well_cells + perf_index_end,
well_cells_.begin() );
well_index_.resize(number_of_perforations_);
std::copy(wells->WI + perf_index_begin,
wells->WI + perf_index_end,
well_index_.begin() );
saturation_table_number_.resize(number_of_perforations_);
std::copy(wells->sat_table_id + perf_index_begin,
wells->sat_table_id + perf_index_end,
saturation_table_number_.begin() );
int perf = 0;
for (const auto& pd : perf_data) {
well_cells_[perf] = pd.cell_index;
well_index_[perf] = pd.connection_transmissibility_factor;
saturation_table_number_[perf] = pd.satnum_id;
++perf;
}
int all_perf = 0;
originalConnectionIndex_.reserve(perf_data.size());
for (const auto connection : well.getConnections()) {
if (connection.state() == Connection::State::OPEN) {
originalConnectionIndex_.push_back(all_perf);
}
++all_perf;
}
assert(originalConnectionIndex_.size() == perf_data.size());
}
// initialization of the completions mapping
@ -205,11 +188,22 @@ namespace Opm
template<typename TypeTag>
WellType
bool
WellInterface<TypeTag>::
wellType() const
isInjector() const
{
return well_type_;
return well_ecl_.isInjector();
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
isProducer() const
{
return well_ecl_.isProducer();
}
@ -781,7 +775,7 @@ namespace Opm
{
// currently, we only updateWellTestState for producers
if (wellType() != PRODUCER) {
if (this->isInjector()) {
return;
}
@ -930,7 +924,8 @@ namespace Opm
bool allCompletionsClosed = true;
const auto& connections = well_ecl_.getConnections();
for (const auto& connection : connections) {
if (!well_test_state.hasCompletion(name(), connection.complnum())) {
if (connection.state() == Connection::State::OPEN
&& !well_test_state.hasCompletion(name(), connection.complnum())) {
allCompletionsClosed = false;
}
}
@ -1166,7 +1161,7 @@ namespace Opm
// we need to get the table number through the parser, in case THP constraint/target is not there.
// When THP control/limit is not active, if available VFP table is provided, we will still need to
// update THP value. However, it will only used for output purpose.
if (well_type_ == PRODUCER) { // producer
if (isProducer()) { // producer
const int table_id = well_ecl_.vfp_table_number();
if (table_id <= 0) {
return false;
@ -1264,10 +1259,12 @@ namespace Opm
const auto& connections = well_ecl_.getConnections();
int perfIdx = 0;
for (const auto& connection : connections) {
if (wellTestState.hasCompletion(name(), connection.complnum())) {
well_index_[perfIdx] = 0.0;
if (connection.state() == Connection::State::OPEN) {
if (wellTestState.hasCompletion(name(), connection.complnum())) {
well_index_[perfIdx] = 0.0;
}
perfIdx++;
}
perfIdx++;
}
}
@ -1295,7 +1292,7 @@ namespace Opm
void
WellInterface<TypeTag>::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger)
{
const auto& connection = well_ecl_.getConnections()[perfIdx];
const auto& connection = well_ecl_.getConnections()[originalConnectionIndex_[perfIdx]];
if (well_ecl_.getDrainageRadius() < 0) {
if (new_well && perfIdx == 0) {
deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero");

316
opm/simulators/wells/WellState.hpp
View File

@ -24,6 +24,8 @@
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/output/data/Wells.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <array>
#include <map>
@ -42,139 +44,54 @@ namespace Opm
typedef std::array< int, 3 > mapentry_t;
typedef std::map< std::string, mapentry_t > WellMapType;
template <class State>
void init(const Wells* wells, const State& state)
{
init(wells, state.pressure());
}
/// Allocate and initialize if wells is non-null.
/// Also tries to give useful initial values to the bhp() and
/// wellRates() fields, depending on controls. The
/// perfRates() field is filled with zero, and perfPress()
/// with -1e100.
void init(const Wells* wells, const std::vector<double>& cellPressures)
void init(const std::vector<double>& cellPressures,
const std::vector<Well>& wells_ecl,
const PhaseUsage& pu,
const std::vector<std::vector<PerforationData>>& well_perf_data,
const SummaryState& summary_state)
{
// clear old name mapping
wellMap_.clear();
wells_.reset( clone_wells( wells ) );
if (wells) {
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
bhp_.resize(nw);
thp_.resize(nw);
well_perf_data_ = well_perf_data;
{
// const int nw = wells->number_of_wells;
const int nw = wells_ecl.size();
// const int np = wells->number_of_phases;
const int np = pu.num_phases;
bhp_.resize(nw, 0.0);
thp_.resize(nw, 0.0);
temperature_.resize(nw, 273.15 + 20); // standard temperature for now
wellrates_.resize(nw * np, 0.0);
int connpos = 0;
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
const Well& well = wells_ecl[w];
// setup wellname -> well index mapping
{
assert( wells->name[ w ] );
std::string name( wells->name[ w ] );
assert( name.size() > 0 );
mapentry_t& wellMapEntry = wellMap_[name];
wellMapEntry[ 0 ] = w;
wellMapEntry[ 1 ] = wells->well_connpos[w];
// also store the number of perforations in this well
wellMapEntry[ 2 ] = num_perf_this_well;
}
// Initialize bhp(), thp(), wellRates().
initSingleWell(cellPressures, w, well, pu, summary_state);
if ( num_perf_this_well == 0 )
{
// No perforations of the well. Initialize to zero.
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = 0.0;
}
bhp_[w] = 0.;
thp_[w] = 0.;
continue;
}
if (well_controls_well_is_stopped(ctrl)) {
// Stopped well:
// 1. Rates: assign zero well rates.
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = 0.0;
}
// 2. Bhp: assign bhp equal to bhp control, if
// applicable, otherwise assign equal to
// first perforation cell pressure.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp_[w] = well_controls_get_current_target( ctrl );
} else {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = cellPressures[first_cell];
}
} else if (well_controls_get_current(ctrl) == -1) {
// Well under group control.
// 1. Rates: assign zero well rates.
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = 0.0;
}
// 2. Bhp: initialize bhp to be a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*cellPressures[first_cell];
} else {
// Open well:
// 1. Rates: initialize well rates to match controls
// if type is SURFACE_RATE. Otherwise, we
// cannot set the correct value here, so we
// assign a small rate with the correct
// sign so that any logic depending on that
// sign will work as expected.
if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
const double rate_target = well_controls_get_current_target(ctrl);
const double * distr = well_controls_get_current_distr( ctrl );
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = rate_target * distr[p];
}
} else {
const double small_rate = 0.0; //1e-14;
const double sign = (wells->type[w] == INJECTOR) ? 1.0 : -1.0;
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = small_rate * sign;
}
}
// 2. Bhp: initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp_[w] = well_controls_get_current_target( ctrl );
} else {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*cellPressures[first_cell];
}
}
// 3. Thp: assign thp equal to thp target/limit, if applicable,
// otherwise keep it zero. Basically, the value should not be used
// in the simulation at all.
const int nwc = well_controls_get_num(ctrl);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
thp_[w] = well_controls_iget_target(ctrl, ctrl_index);
break;
}
}
// Setup wellname -> well index mapping.
const int num_perf_this_well = well_perf_data[w].size();
std::string name = well.name();
assert( name.size() > 0 );
mapentry_t& wellMapEntry = wellMap_[name];
wellMapEntry[ 0 ] = w;
wellMapEntry[ 1 ] = connpos;
wellMapEntry[ 2 ] = num_perf_this_well;
connpos += num_perf_this_well;
}
// The perforation rates and perforation pressures are
// not expected to be consistent with bhp_ and wellrates_
// after init().
perfrates_.resize(wells->well_connpos[nw], 0.0);
perfpress_.resize(wells->well_connpos[nw], -1e100);
perfrates_.resize(connpos, 0.0);
perfpress_.resize(connpos, -1e100);
}
}
@ -264,14 +181,11 @@ namespace Opm
well.rates.set( rt::gas, wv[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ] );
}
const int num_perf_well = this->wells_->well_connpos[ well_index + 1 ]
- this->wells_->well_connpos[ well_index ];
const int num_perf_well = this->well_perf_data_[well_index].size();
well.connections.resize(num_perf_well);
for( int i = 0; i < num_perf_well; ++i ) {
const auto wi = this->wells_->well_connpos[ well_index ] + i;
const auto active_index = this->wells_->well_cells[ wi ];
const auto active_index = this->well_perf_data_[well_index][i].cell_index;
auto& connection = well.connections[ i ];
connection.index = globalCellIdxMap[active_index];
connection.pressure = this->perfPress()[ itr.second[1] + i ];
@ -284,32 +198,10 @@ namespace Opm
}
virtual ~WellState() {}
virtual ~WellState() = default;
WellState() = default;
WellState( const WellState& rhs ) :
bhp_( rhs.bhp_ ),
thp_( rhs.thp_ ),
temperature_( rhs.temperature_ ),
wellrates_( rhs.wellrates_ ),
perfrates_( rhs.perfrates_ ),
perfpress_( rhs.perfpress_ ),
wellMap_( rhs.wellMap_ ),
wells_( clone_wells( rhs.wells_.get() ) )
{}
WellState& operator=( const WellState& rhs ) {
this->bhp_ = rhs.bhp_;
this->thp_ = rhs.thp_;
this->temperature_ = rhs.temperature_;
this->wellrates_ = rhs.wellrates_;
this->perfrates_ = rhs.perfrates_;
this->perfpress_ = rhs.perfpress_;
this->wellMap_ = rhs.wellMap_;
this->wells_.reset( clone_wells( rhs.wells_.get() ) );
return *this;
}
WellState(const WellState& rhs) = default;
WellState& operator=(const WellState& rhs) = default;
private:
std::vector<double> bhp_;
@ -321,11 +213,139 @@ namespace Opm
WellMapType wellMap_;
void initSingleWell(const std::vector<double>& cellPressures,
const int w,
const Well& well,
const PhaseUsage& pu,
const SummaryState& summary_state)
{
assert(well.isInjector() || well.isProducer());
// Set default zero initial well rates.
// May be overwritten below.
const int np = pu.num_phases;
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = 0.0;
}
const int num_perf_this_well = well_perf_data_[w].size();
if ( num_perf_this_well == 0 ) {
// No perforations of the well. Initialize to zero.
bhp_[w] = 0.;
thp_[w] = 0.;
return;
}
const auto inj_controls = well.isInjector() ? well.injectionControls(summary_state) : Well::InjectionControls(0);
const auto prod_controls = well.isProducer() ? well.productionControls(summary_state) : Well::ProductionControls(0);
const bool is_bhp = well.isInjector() ? (inj_controls.cmode == Well::InjectorCMode::BHP)
: (prod_controls.cmode == Well::ProducerCMode::BHP);
const double bhp_limit = well.isInjector() ? inj_controls.bhp_limit : prod_controls.bhp_limit;
const bool is_grup = well.isInjector() ? (inj_controls.cmode == Well::InjectorCMode::GRUP)
: (prod_controls.cmode == Well::ProducerCMode::GRUP);
const double inj_surf_rate = well.isInjector() ? inj_controls.surface_rate : 0.0; // To avoid a "maybe-uninitialized" warning.
if (well.getStatus() == Well::Status::STOP) {
// Stopped well:
// 1. Rates: zero well rates.
// 2. Bhp: assign bhp equal to bhp control, if
// applicable, otherwise assign equal to
// first perforation cell pressure.
if (is_bhp) {
bhp_[w] = bhp_limit;
} else {
const int first_cell = well_perf_data_[w][0].cell_index;
bhp_[w] = cellPressures[first_cell];
}
} else if (is_grup) {
// Well under group control.
// 1. Rates: zero well rates.
// 2. Bhp: initialize bhp to be a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
const int first_cell = well_perf_data_[w][0].cell_index;
const double safety_factor = well.isInjector() ? 1.01 : 0.99;
bhp_[w] = safety_factor*cellPressures[first_cell];
} else {
// Open well, under own control:
// 1. Rates: initialize well rates to match
// controls if type is ORAT/GRAT/WRAT
// (producer) or RATE (injector).
// Otherwise, we cannot set the correct
// value here and initialize to zero rate.
if (well.isInjector()) {
if (inj_controls.cmode == Well::InjectorCMode::RATE) {
switch (inj_controls.injector_type) {
case Well::InjectorType::WATER:
assert(pu.phase_used[BlackoilPhases::Aqua]);
wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
break;
case Well::InjectorType::GAS:
assert(pu.phase_used[BlackoilPhases::Vapour]);
wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
break;
case Well::InjectorType::OIL:
assert(pu.phase_used[BlackoilPhases::Liquid]);
wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
break;
case Well::InjectorType::MULTI:
// Not currently handled, keep zero init.
break;
}
} else {
// Keep zero init.
}
} else {
assert(well.isProducer());
// Note negative rates for producing wells.
switch (prod_controls.cmode) {
case Well::ProducerCMode::ORAT:
assert(pu.phase_used[BlackoilPhases::Liquid]);
wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
break;
case Well::ProducerCMode::WRAT:
assert(pu.phase_used[BlackoilPhases::Aqua]);
wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
break;
case Well::ProducerCMode::GRAT:
assert(pu.phase_used[BlackoilPhases::Vapour]);
wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
break;
default:
// Keep zero init.
break;
}
}
// 2. Bhp: initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
if (is_bhp) {
bhp_[w] = bhp_limit;
} else {
const int first_cell = well_perf_data_[w][0].cell_index;
const double safety_factor = well.isInjector() ? 1.01 : 0.99;
bhp_[w] = safety_factor*cellPressures[first_cell];
}
}
// 3. Thp: assign thp equal to thp target/limit, if such a limit exists,
// otherwise keep it zero.
const bool has_thp = well.isInjector() ? inj_controls.hasControl(Well::InjectorCMode::THP)
: prod_controls.hasControl(Well::ProducerCMode::THP);
const double thp_limit = well.isInjector() ? inj_controls.thp_limit : prod_controls.thp_limit;
if (has_thp) {
thp_[w] = thp_limit;
}
}
protected:
struct wdel {
void operator()( Wells* w ) { destroy_wells( w ); }
};
std::unique_ptr< Wells, wdel > wells_;
std::vector<std::vector<PerforationData>> well_perf_data_;
};
} // namespace Opm

177
opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
View File

@ -52,6 +52,8 @@ namespace Opm
public:
typedef BaseType :: WellMapType WellMapType;
virtual ~WellStateFullyImplicitBlackoil() = default;
// TODO: same definition with WellInterface, eventually they should go to a common header file.
static const int Water = BlackoilPhases::Aqua;
static const int Oil = BlackoilPhases::Liquid;
@ -67,29 +69,29 @@ namespace Opm
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls
void init(const Wells* wells,
const std::vector<double>& cellPressures,
void init(const std::vector<double>& cellPressures,
const Schedule& schedule,
const std::vector<Well>& wells_ecl,
const int report_step,
const WellStateFullyImplicitBlackoil* prevState,
const PhaseUsage& pu)
const PhaseUsage& pu,
const std::vector<std::vector<PerforationData>>& well_perf_data,
const SummaryState& summary_state)
{
// call init on base class
BaseType :: init(wells, cellPressures);
BaseType :: init(cellPressures, wells_ecl, pu, well_perf_data, summary_state);
// if there are no well, do nothing in init
if (wells == 0) {
return;
}
const int nw = wells->number_of_wells;
const int nw = wells_ecl.size();
if( nw == 0 ) return ;
// Initialize perfphaserates_, which must be done here.
const int np = wells->number_of_phases;
const int nperf = wells->well_connpos[nw];
const int np = pu.num_phases;
int nperf = 0;
for (const auto& wpd : well_perf_data) {
nperf += wpd.size();
}
well_reservoir_rates_.resize(nw * np, 0.0);
well_dissolved_gas_rates_.resize(nw, 0.0);
@ -107,23 +109,9 @@ namespace Opm
const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+ ScheduleEvents::PRODUCTION_UPDATE
+ ScheduleEvents::INJECTION_UPDATE;
for (int w = 0; w <nw; ++w) {
const int nw_wells_ecl = wells_ecl.size();
int index_well_ecl = 0;
const std::string well_name(wells->name[w]);
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_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
}
effective_events_occurred_[w] = (schedule.hasWellEvent(well_name, effective_events_mask, report_step) );
for (int w = 0; w < nw; ++w) {
effective_events_occurred_[w]
= schedule.hasWellEvent(wells_ecl[w].name(), effective_events_mask, report_step);
}
} // end of if (!well_ecl.empty() )
@ -139,21 +127,20 @@ namespace Opm
perf_skin_pressure_.clear();
perf_skin_pressure_.resize(nperf, 0.0);
int connpos = 0;
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
// Open well: Initialize perfphaserates_ to well
// Initialize perfphaserates_ to well
// rates divided by the number of perforations.
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
if (well_controls_well_is_open(ctrl)) {
const int num_perf_this_well = well_perf_data[w].size();
for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) {
if (wells_ecl[w].getStatus() == Well::Status::OPEN) {
for (int p = 0; p < np; ++p) {
perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
}
}
perfPress()[perf] = cellPressures[wells->well_cells[perf]];
perfPress()[perf] = cellPressures[well_perf_data[w][perf-connpos].cell_index];
}
connpos += num_perf_this_well;
}
current_injection_controls_.resize(nw);
@ -166,13 +153,14 @@ namespace Opm
// intialize wells that have been there before
// order may change so the mapping is based on the well name
if(prevState && !prevState->wellMap().empty()) {
typedef typename WellMapType :: const_iterator const_iterator;
const_iterator end = prevState->wellMap().end();
if (prevState && !prevState->wellMap().empty()) {
connpos = 0;
auto end = prevState->wellMap().end();
for (int w = 0; w < nw; ++w) {
const std::string name( wells->name[ w ] );
const_iterator it = prevState->wellMap().find( name );
if( it != end )
const Well& well = wells_ecl[w];
const int num_perf_this_well = well_perf_data[w].size();
auto it = prevState->wellMap().find(well.name());
if ( it != end )
{
const int oldIndex = (*it).second[ 0 ];
const int newIndex = w;
@ -210,19 +198,18 @@ namespace Opm
// perfPhaseRates
const int oldPerf_idx_beg = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ];
const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
// copy perforation rates when the number of perforations is equal,
// otherwise initialize perfphaserates to well rates divided by the number of perforations.
if( num_perf_old_well == num_perf_this_well )
{
int old_perf_phase_idx = oldPerf_idx_beg *np;
for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
for (int perf_phase_idx = connpos*np;
perf_phase_idx < (connpos + num_perf_this_well)*np; ++perf_phase_idx, ++old_perf_phase_idx )
{
perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
}
} else {
for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) {
for (int p = 0; p < np; ++p) {
perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well);
}
@ -232,8 +219,7 @@ namespace Opm
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
{
perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ];
}
@ -243,8 +229,7 @@ namespace Opm
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
{
perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ];
}
@ -258,8 +243,7 @@ namespace Opm
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
{
perf_water_throughput_[ perf ] = prevState->perfThroughput()[ oldPerf_idx ];
perf_skin_pressure_[ perf ] = prevState->perfSkinPressure()[ oldPerf_idx ];
@ -269,21 +253,16 @@ namespace Opm
}
}
// If in the new step, there is no THP related target/limit anymore, its thp value should be
// set to zero.
const WellControls* ctrl = wells->ctrls[w];
const int nwc = well_controls_get_num(ctrl);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
break;
}
}
// not finding any thp related control/limits
if (ctrl_index == nwc) {
thp()[w] = 0.;
const bool has_thp = well.isInjector() ? well.injectionControls(summary_state).hasControl(Well::InjectorCMode::THP)
: well.productionControls(summary_state).hasControl(Well::ProducerCMode::THP);
if (!has_thp) {
thp()[w] = 0.0;
}
// Increment connection position offset.
connpos += num_perf_this_well;
}
}
@ -303,14 +282,19 @@ namespace Opm
}
void resize(const Wells* wells, const std::vector<Well>& wells_ecl, const Schedule& schedule,
const bool handle_ms_well, const size_t numCells, const PhaseUsage& pu)
void resize(const std::vector<Well>& wells_ecl,
const Schedule& schedule,
const bool handle_ms_well,
const size_t numCells,
const PhaseUsage& pu,
const std::vector<std::vector<PerforationData>>& well_perf_data,
const SummaryState& summary_state)
{
const std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
init(wells, tmp, schedule, wells_ecl, 0, nullptr, pu);
init(tmp, schedule, wells_ecl, 0, nullptr, pu, well_perf_data, summary_state);
if (handle_ms_well) {
initWellStateMSWell(wells, wells_ecl, pu, nullptr);
initWellStateMSWell(wells_ecl, pu, nullptr);
}
}
@ -434,13 +418,6 @@ namespace Opm
phs.at( pu.phase_pos[Gas] ) = rt::gas;
}
if (pu.has_solvent) {
// add solvent component
for( int w = 0; w < nw; ++w ) {
res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) );
}
}
/* this is a reference or example on **how** to convert from
* WellState to something understood by opm-output. it is intended
* to be properly implemented and maintained as a part of
@ -491,10 +468,14 @@ namespace Opm
well.rates.set( rt::well_potential_gas, this->well_potentials_[well_rate_index + pu.phase_pos[Gas]] );
}
if ( pu.has_solvent ) {
well.rates.set( rt::solvent, solventWellRate(w) );
}
well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] );
well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] );
int local_comp_index = 0;
size_t local_comp_index = 0;
for( auto& comp : well.connections) {
const auto rates = this->perfPhaseRates().begin()
+ (np * wt.second[ 1 ])
@ -505,7 +486,7 @@ namespace Opm
comp.rates.set( phs[ i ], *(rates + i) );
}
}
assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]);
assert(local_comp_index == this->well_perf_data_[w].size());
const auto nseg = this->numSegments(w);
for (auto seg_ix = 0*nseg; seg_ix < nseg; ++seg_ix) {
@ -520,11 +501,11 @@ namespace Opm
/// init the MS well related.
void initWellStateMSWell(const Wells* wells, const std::vector<Well>& wells_ecl,
void initWellStateMSWell(const std::vector<Well>& wells_ecl,
const PhaseUsage& pu, const WellStateFullyImplicitBlackoil* prev_well_state)
{
// still using the order in wells
const int nw = wells->number_of_wells;
const int nw = wells_ecl.size();
if (nw == 0) {
return;
}
@ -538,24 +519,12 @@ namespace Opm
seg_number_.clear();
nseg_ = 0;
int connpos = 0;
// in the init function, the well rates and perforation rates have been initialized or copied from prevState
// what we do here, is to set the segment rates and perforation rates
for (int w = 0; w < nw; ++w) {
const int nw_wells_ecl = wells_ecl.size();
int index_well_ecl = 0;
const std::string well_name(wells->name[w]);
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_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
}
const auto& well_ecl = wells_ecl[index_well_ecl];
const auto& well_ecl = wells_ecl[w];
int num_perf_this_well = well_perf_data_[w].size();
top_segment_index_.push_back(nseg_);
if ( !well_ecl.isMultiSegment() ) { // not multi-segment well
nseg_ += 1;
@ -571,7 +540,6 @@ namespace Opm
const WellConnections& completion_set = well_ecl.getConnections();
// number of segment for this single well
const int well_nseg = segment_set.size();
// const int nperf = completion_set.size();
int n_activeperf = 0;
nseg_ += well_nseg;
for (auto segID = 0*well_nseg; segID < well_nseg; ++segID) {
@ -605,8 +573,8 @@ namespace Opm
// for the segrates_, now it becomes a recursive solution procedure.
{
const int np = numPhases();
const int start_perf = wells->well_connpos[w];
const int start_perf_next_well = wells->well_connpos[w + 1];
const int start_perf = connpos;
const int start_perf_next_well = connpos + num_perf_this_well;
// make sure the information from wells_ecl consistent with wells
assert(n_activeperf == (start_perf_next_well - start_perf));
@ -640,7 +608,7 @@ namespace Opm
// top segment is always the first one, and its pressure is the well bhp
segpress_.push_back(bhp()[w]);
const int top_segment = top_segment_index_[w];
const int start_perf = wells->well_connpos[w];
const int start_perf = connpos;
for (int seg = 1; seg < well_nseg; ++seg) {
if ( !segment_perforations[seg].empty() ) {
const int first_perf = segment_perforations[seg][0];
@ -654,6 +622,7 @@ namespace Opm
}
}
}
connpos += num_perf_this_well;
}
assert(int(segpress_.size()) == nseg_);
assert(int(segrates_.size()) == nseg_ * numPhases() );
@ -663,8 +632,7 @@ namespace Opm
const auto& end = prev_well_state->wellMap().end();
const int np = numPhases();
for (int w = 0; w < nw; ++w) {
const std::string name( wells->name[w] );
const auto& it = prev_well_state->wellMap().find( name );
const auto& it = prev_well_state->wellMap().find( wells_ecl[w].name() );
if (it != end) { // the well is found in the prev_well_state
// TODO: the well with same name can change a lot, like they might not have same number of segments
@ -736,8 +704,13 @@ namespace Opm
/// One rate pr well
double solventWellRate(const int w) const {
int connpos = 0;
for (int iw = 0; iw < w; ++iw) {
connpos += this->well_perf_data_[iw].size();
}
double solvent_well_rate = 0.0;
for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
const int endperf = connpos + this->well_perf_data_[w].size();
for (int perf = connpos; perf < endperf; ++perf ) {
solvent_well_rate += perfRateSolvent_[perf];
}
return solvent_well_rate;

48
tests/test_stoppedwells.cpp
View File

@ -27,12 +27,6 @@
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/grid/GridManager.hpp>
using namespace Opm;
@ -42,54 +36,24 @@ BOOST_AUTO_TEST_CASE(TestStoppedWells)
Opm::Parser parser;
Opm::Deck deck(parser.parseFile(filename));
Opm::EclipseState eclipseState(deck);
Opm::GridManager vanguard(eclipseState.getInputGrid());
const auto& grid = eclipseState.getInputGrid();
const TableManager table ( deck );
const Eclipse3DProperties eclipseProperties ( deck , table, grid);
const Opm::Runspec runspec (deck);
const Schedule sched(deck, grid, eclipseProperties, runspec);
Opm::SummaryState summaryState(std::chrono::system_clock::from_time_t(sched.getStartTime()));
double target_surfacerate_inj;
double target_surfacerate_prod;
const std::vector<double> pressure = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
// Both wells are open in the first schedule step
{
Opm::WellsManager wellsManager(eclipseState , sched, summaryState, 0, *vanguard.c_grid());
const Wells* wells = wellsManager.c_wells();
const struct WellControls* ctrls0 = wells->ctrls[0];
const struct WellControls* ctrls1 = wells->ctrls[1];
BOOST_CHECK(well_controls_well_is_open(ctrls0));
BOOST_CHECK(well_controls_well_is_open(ctrls1));
target_surfacerate_inj = well_controls_iget_target(ctrls0 , 0);
target_surfacerate_prod = well_controls_iget_target(ctrls1 , 0);
WellState wellstate;
wellstate.init(wells, pressure);
const std::vector<double> wellrates = wellstate.wellRates();
BOOST_CHECK_EQUAL (target_surfacerate_inj, wellrates[2]); // Gas injector
BOOST_CHECK_EQUAL (target_surfacerate_prod, wellrates[4]); // Oil target rate
auto wells = sched.getWells(0);
BOOST_CHECK(wells[0].getStatus() == Opm::Well::Status::OPEN);
BOOST_CHECK(wells[1].getStatus() == Opm::Well::Status::OPEN);
}
// The injector is stopped
{
Opm::WellsManager wellsManager(eclipseState, sched, summaryState, 1 , *vanguard.c_grid());
const Wells* wells = wellsManager.c_wells();
const struct WellControls* ctrls0 = wells->ctrls[0];
const struct WellControls* ctrls1 = wells->ctrls[1];
BOOST_CHECK(well_controls_well_is_stopped(ctrls0)); // injector is stopped
BOOST_CHECK(well_controls_well_is_open(ctrls1));
WellState wellstate;
wellstate.init(wells, pressure);
const std::vector<double> wellrates = wellstate.wellRates();
BOOST_CHECK_EQUAL (0, wellrates[2]); // Gas injector
BOOST_CHECK_EQUAL (target_surfacerate_prod, wellrates[4]); // Oil target rate
auto wells = sched.getWells(1);
BOOST_CHECK(wells[0].getStatus() == Opm::Well::Status::STOP);
BOOST_CHECK(wells[1].getStatus() == Opm::Well::Status::OPEN);
}
}

58
tests/test_wellmodel.cpp
View File

@ -77,36 +77,13 @@ struct SetupTest {
schedule.reset( new Opm::Schedule(deck, ecl_state->getInputGrid(), eclipseProperties, runspec));
summaryState.reset( new Opm::SummaryState(std::chrono::system_clock::from_time_t(schedule->getStartTime())));
}
// Create grid.
const std::vector<double>& porv =
ecl_state->get3DProperties().getDoubleGridProperty("PORV").getData();
Opm::GridManager gm(ecl_state->getInputGrid(), porv);
const Grid& grid = *(gm.c_grid());
current_timestep = 0;
// Create wells.
wells_manager.reset(new Opm::WellsManager(*ecl_state,
*schedule,
*summaryState,
current_timestep,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::dimensions(grid),
Opm::UgGridHelpers::cell2Faces(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
false,
std::unordered_set<std::string>() ) );
};
std::unique_ptr<const Opm::WellsManager> wells_manager;
std::unique_ptr<const Opm::EclipseState> ecl_state;
std::unique_ptr<const Opm::Schedule> schedule;
std::unique_ptr<Opm::SummaryState> summaryState;
std::vector<std::vector<Opm::PerforationData>> well_perf_data;
int current_timestep;
};
@ -132,7 +109,6 @@ BOOST_GLOBAL_FIXTURE(GlobalFixture);
BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const SetupTest setup_test;
const Wells* wells = setup_test.wells_manager->c_wells();
const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep);
BOOST_CHECK_EQUAL( wells_ecl.size(), 2);
const Opm::Well& well = wells_ecl[1];
@ -149,36 +125,24 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
rateConverter.reset(new RateConverterType (phaseUsage,
std::vector<int>(10, 0)));
const int pvtIdx = 0;
const int num_comp = wells->number_of_phases;
Opm::PerforationData dummy;
std::vector<Opm::PerforationData> pdata(well.getConnections().size(), dummy);
BOOST_CHECK_THROW( StandardWell( well, -1, wells, param, *rateConverter, pvtIdx, num_comp), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr , param, *rateConverter, pvtIdx, num_comp), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, -1, param, *rateConverter, 0, 3, 3, 0, 0, pdata), std::invalid_argument);
}
BOOST_AUTO_TEST_CASE(TestBehavoir) {
const SetupTest setup_test;
const Wells* wells_struct = setup_test.wells_manager->c_wells();
const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep);
const int current_timestep = setup_test.current_timestep;
std::vector<std::unique_ptr<const StandardWell> > wells;
{
const int nw = wells_struct ? (wells_struct->number_of_wells) : 0;
const int nw = wells_ecl.size();
const Opm::BlackoilModelParametersEbos<TTAG(EclFlowProblem)> param;
for (int w = 0; w < nw; ++w) {
const std::string well_name(wells_struct->name[w]);
size_t index_well = 0;
for (; index_well < wells_ecl.size(); ++index_well) {
if (well_name == wells_ecl[index_well].name()) {
break;
}
}
// we should always be able to find the well in wells_ecl
BOOST_CHECK(index_well != wells_ecl.size());
// For the conversion between the surface volume rate and resrevoir voidage rate
typedef Opm::BlackOilFluidSystem<double> FluidSystem;
using RateConverterType = Opm::RateConverter::
@ -191,11 +155,9 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
// Compute reservoir volumes for RESV controls.
rateConverter.reset(new RateConverterType (phaseUsage,
std::vector<int>(10, 0)));
const int pvtIdx = 0;
const int num_comp = wells_struct->number_of_phases;
wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param, *rateConverter, pvtIdx, num_comp) );
Opm::PerforationData dummy;
std::vector<Opm::PerforationData> pdata(wells_ecl[w].getConnections().size(), dummy);
wells.emplace_back(new StandardWell(wells_ecl[w], current_timestep, param, *rateConverter, 0, 3, 3, w, 0, pdata) );
}
}
@ -203,7 +165,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{
const auto& well = wells[0];
BOOST_CHECK_EQUAL(well->name(), "PROD1");
BOOST_CHECK(well->wellType() == PRODUCER);
BOOST_CHECK(well->isProducer());
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4);
}
@ -212,7 +174,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{
const auto& well = wells[1];
BOOST_CHECK_EQUAL(well->name(), "INJE1");
BOOST_CHECK(well->wellType() == INJECTOR);
BOOST_CHECK(well->isInjector());
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4);
}

58
tests/test_wellstatefullyimplicitblackoil.cpp
View File

@ -53,13 +53,60 @@ struct Setup
, grid (es.getInputGrid())
, sched(deck, es)
, st(std::chrono::system_clock::from_time_t(sched.getStartTime()))
{}
{
initWellPerfData();
}
void initWellPerfData()
{
const auto& wells = sched.getWells(0);
const auto& cartDims = Opm::UgGridHelpers::cartDims(*grid.c_grid());
const int* compressed_to_cartesian = Opm::UgGridHelpers::globalCell(*grid.c_grid());
std::vector<int> cartesian_to_compressed(cartDims[0] * cartDims[1] * cartDims[2], -1);
for (int ii = 0; ii < Opm::UgGridHelpers::numCells(*grid.c_grid()); ++ii) {
cartesian_to_compressed[compressed_to_cartesian[ii]] = ii;
}
well_perf_data.resize(wells.size());
int well_index = 0;
for (const auto& well : wells) {
well_perf_data[well_index].clear();
well_perf_data[well_index].reserve(well.getConnections().size());
for (const auto& completion : well.getConnections()) {
if (completion.state() == Opm::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 {
Opm::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() != Opm::Connection::State::SHUT) {
OPM_THROW(std::runtime_error,
"Completion state: " << Opm::Connection::State2String(completion.state()) << " not handled");
}
}
}
++well_index;
}
}
Opm::EclipseState es;
Opm::PhaseUsage pu;
Opm::GridManager grid;
Opm::Schedule sched;
Opm::SummaryState st;
std::vector<std::vector<Opm::PerforationData>> well_perf_data;
};
namespace {
@ -72,14 +119,11 @@ namespace {
std::vector<double>(setup.grid.c_grid()->number_of_cells,
100.0*Opm::unit::barsa);
const Opm::WellsManager wmgr{setup.es, setup.sched, setup.st, timeStep, *setup.grid.c_grid()};
state.init(wmgr.c_wells(), cpress, setup.sched,
state.init(cpress, setup.sched,
setup.sched.getWells(timeStep),
timeStep, nullptr, setup.pu);
timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st);
state.initWellStateMSWell(wmgr.c_wells(),
setup.sched.getWells(timeStep),
state.initWellStateMSWell(setup.sched.getWells(timeStep),
setup.pu, nullptr);
return state;