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

Merge pull request #2164 from atgeirr/remove-wells-struct-rebased

Browse Source 
Avoid using the Wells struct.
This commit is contained in:
Bård Skaflestad 2019-11-25 09:56:06 +01:00 committed by GitHub
commit 334acd18ad
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
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/gatherDeferredLogger.hpp
opm/simulators/utils/moduleVersion.hpp opm/simulators/utils/moduleVersion.hpp
opm/simulators/utils/ParallelRestart.hpp opm/simulators/utils/ParallelRestart.hpp
opm/simulators/wells/PerforationData.hpp
opm/simulators/wells/RateConverter.hpp opm/simulators/wells/RateConverter.hpp
opm/simulators/wells/SimFIBODetails.hpp opm/simulators/wells/SimFIBODetails.hpp
opm/simulators/wells/WellConnectionAuxiliaryModule.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.h>
#include <opm/core/wells/WellCollection.hpp> #include <opm/core/wells/WellCollection.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp> #include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <opm/simulators/wells/VFPInjProperties.hpp> #include <opm/simulators/wells/VFPInjProperties.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp> #include <opm/simulators/wells/VFPProdProperties.hpp>
#include <opm/simulators/flow/countGlobalCells.hpp> #include <opm/simulators/flow/countGlobalCells.hpp>
@ -233,8 +234,10 @@ namespace Opm {
protected: protected:
Simulator& ebosSimulator_; Simulator& ebosSimulator_;
std::unique_ptr<WellsManager> wells_manager_;
std::vector< Well > wells_ecl_; std::vector< Well > wells_ecl_;
std::vector< std::vector<PerforationData> > well_perf_data_;
std::vector<int> first_perf_index_;
bool wells_active_; bool wells_active_;
@ -246,8 +249,10 @@ namespace Opm {
std::vector<bool> is_cell_perforated_; std::vector<bool> is_cell_perforated_;
void initializeWellPerfData();
// create the well container // 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; 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 // used to better efficiency of calcuation
mutable BVector scaleAddRes_; mutable BVector scaleAddRes_;
const Wells* wells() const { return wells_manager_->c_wells(); }
const Grid& grid() const const Grid& grid() const
{ return ebosSimulator_.vanguard().grid(); } { return ebosSimulator_.vanguard().grid(); }
@ -373,7 +376,7 @@ namespace Opm {
WellStateFullyImplicitBlackoil& state ) const; WellStateFullyImplicitBlackoil& state ) const;
// whether there exists any multisegment well open on this process // 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; const Well& getWellEcl(const std::string& well_name) const;

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

@ -37,6 +37,19 @@ namespace Opm {
// Create the guide rate container. // Create the guide rate container.
guideRate_.reset(new GuideRate (ebosSimulator_.vanguard().schedule())); 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> template<typename TypeTag>
@ -46,31 +59,15 @@ namespace Opm {
{ {
const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState(); const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState();
gravity_ = ebosSimulator_.problem().gravity()[2];
extractLegacyCellPvtRegionIndex_(); extractLegacyCellPvtRegionIndex_();
extractLegacyDepth_(); extractLegacyDepth_();
phase_usage_ = phaseUsageFromDeck(eclState); 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]; gravity_ = ebosSimulator_.problem().gravity()[2];
extractLegacyCellPvtRegionIndex_();
extractLegacyDepth_();
initial_step_ = true; 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 // add the eWoms auxiliary module for the wells to the list
ebosSimulator_.model().addAuxiliaryModule(this); ebosSimulator_.model().addAuxiliaryModule(this);
@ -210,24 +207,18 @@ namespace Opm {
Opm::DeferredLogger local_deferredLogger; Opm::DeferredLogger local_deferredLogger;
const Grid& grid = ebosSimulator_.vanguard().grid(); 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(); const auto& summaryState = ebosSimulator_.vanguard().summaryState();
wells_ecl_ = schedule().getWells(timeStepIdx); // Make wells_ecl_ contain only this partition's non-shut wells.
{
// Create wells and well state. const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
wells_manager_.reset( new WellsManager (eclState, auto is_shut_or_defunct = [&defunct_well_names](const Well& well) {
schedule(), return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end());
summaryState, };
timeStepIdx, auto w = schedule().getWells(timeStepIdx);
Opm::UgGridHelpers::numCells(grid), w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end());
Opm::UgGridHelpers::globalCell(grid), wells_ecl_.swap(w);
Opm::UgGridHelpers::cartDims(grid), }
Opm::UgGridHelpers::dimensions(grid), initializeWellPerfData();
Opm::UgGridHelpers::cell2Faces(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
grid.comm().size() > 1,
defunct_well_names) );
// Wells are active if they are active wells on at least // Wells are active if they are active wells on at least
// one process. // one process.
@ -269,35 +260,22 @@ namespace Opm {
} }
cellPressures[cellIdx] = perf_pressure; 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 // handling MS well related
if (param_.use_multisegment_well_&& anyMSWellOpenLocal(wells())) { // if we use MultisegmentWell model if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
well_state_.initWellStateMSWell(wells(), wells_ecl_, phase_usage_, &previous_well_state_); 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) { for (int w = 0; w <nw; ++w) {
const int nw_wells_ecl = wells_ecl_.size(); const auto& well = wells_ecl_[w];
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 uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+ ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::PRODUCTION_UPDATE
+ ScheduleEvents::INJECTION_UPDATE + ScheduleEvents::INJECTION_UPDATE
+ ScheduleEvents::NEW_WELL; + ScheduleEvents::NEW_WELL;
if(!schedule().hasWellEvent(well_name, effective_events_mask, timeStepIdx)) if(!schedule().hasWellEvent(well.name(), effective_events_mask, timeStepIdx))
continue; continue;
if (well.isProducer()) { if (well.isProducer()) {
@ -347,7 +325,7 @@ namespace Opm {
wellTesting(reportStepIdx, simulationTime, local_deferredLogger); wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
// create the well container // create the well container
well_container_ = createWellContainer(reportStepIdx, local_deferredLogger); well_container_ = createWellContainer(reportStepIdx);
// do the initialization for all the wells // do the initialization for all the wells
// TODO: to see whether we can postpone of the intialization of the well containers to // TODO: to see whether we can postpone of the intialization of the well containers to
@ -400,9 +378,16 @@ namespace Opm {
//compute well guideRates //compute well guideRates
const int np = numPhases(); const int np = numPhases();
for (const auto& well : well_container_) { for (const auto& well : well_container_) {
const double& oilpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Liquid]]; const auto wpot = well_state_.wellPotentials().data() + well->indexOfWell()*np;
const double& gaspot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Vapour]]; const double oilpot = (phase_usage_.phase_used[BlackoilPhases::Liquid] > 0)
const double& waterpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Aqua]]; ? 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); guideRate_->compute(well->name(), reportStepIdx, simulationTime, oilpot, gaspot, waterpot);
} }
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
@ -470,7 +455,7 @@ namespace Opm {
Opm::DeferredLogger local_deferredLogger; Opm::DeferredLogger local_deferredLogger;
for (const auto& well : well_container_) { 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_); well->updateWaterThroughput(dt, well_state_);
} }
} }
@ -538,8 +523,6 @@ namespace Opm {
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
initFromRestartFile(const RestartValue& restartValues) initFromRestartFile(const RestartValue& restartValues)
{ {
const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames();
// The restart step value is used to identify wells present at the given // 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 // 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 // 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 int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0);
const auto& summaryState = ebosSimulator_.vanguard().summaryState(); const auto& summaryState = ebosSimulator_.vanguard().summaryState();
WellsManager wellsmanager(eclState(), // Make wells_ecl_ contain only this partition's non-shut wells.
schedule(), {
summaryState, const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames();
report_step, auto is_shut_or_defunct = [&defunct_well_names](const Well& well) {
Opm::UgGridHelpers::numCells(grid()), return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end());
Opm::UgGridHelpers::globalCell(grid()), };
Opm::UgGridHelpers::cartDims(grid()), auto w = schedule().getWells(report_step);
Opm::UgGridHelpers::dimensions(grid()), w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end());
Opm::UgGridHelpers::cell2Faces(grid()), wells_ecl_.swap(w);
Opm::UgGridHelpers::beginFaceCentroids(grid()), }
grid().comm().size() > 1,
defunctWellNames);
const Wells* wells = wellsmanager.c_wells(); initializeWellPerfData();
wells_ecl_ = schedule().getWells(report_step); const int nw = wells_ecl_.size();
const int nw = wells->number_of_wells;
if (nw > 0) { if (nw > 0) {
const auto phaseUsage = phaseUsageFromDeck(eclState()); const auto phaseUsage = phaseUsageFromDeck(eclState());
const size_t numCells = Opm::UgGridHelpers::numCells(grid()); const size_t numCells = Opm::UgGridHelpers::numCells(grid());
const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal(wells)); const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal());
well_state_.resize(wells, wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage); // Resize for restart step 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_); wellsToState(restartValues.wells, phaseUsage, handle_ms_well, well_state_);
} }
@ -578,20 +557,7 @@ namespace Opm {
const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST(); const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST();
if (true || ecl_compatible_rst) { // always set the control from the schedule if (true || ecl_compatible_rst) { // always set the control from the schedule
for (int w = 0; w <nw; ++w) { for (int w = 0; w <nw; ++w) {
const int nw_wells_ecl = wells_ecl_.size(); const auto& well = wells_ecl_[w];
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];
if (well.isProducer()) { if (well.isProducer()) {
const auto controls = well.productionControls(summaryState); 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> template<typename TypeTag>
std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr > std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger) createWellContainer(const int time_step)
{ {
std::vector<WellInterfacePtr> well_container; std::vector<WellInterfacePtr> well_container;
@ -624,27 +639,9 @@ namespace Opm {
if (nw > 0) { if (nw > 0) {
well_container.reserve(nw); 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) { for (int w = 0; w < nw; ++w) {
const std::string well_name = std::string(wells()->name[w]); const Well& well_ecl = wells_ecl_[w];
const std::string& well_name = well_ecl.name();
// 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];
// A new WCON keywords can re-open a well that was closed/shut due to Physical limit // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
if ( wellTestState_.hasWellClosed(well_name)) { if ( wellTestState_.hasWellClosed(well_name)) {
@ -663,7 +660,6 @@ namespace Opm {
} }
} }
// TODO: should we do this for all kinds of closing reasons? // TODO: should we do this for all kinds of closing reasons?
// something like wellTestState_.hasWell(well_name)? // something like wellTestState_.hasWell(well_name)?
bool wellIsStopped = false; bool wellIsStopped = false;
@ -688,15 +684,31 @@ namespace Opm {
} }
// Use the pvtRegionIdx from the top cell // 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]; const int pvtreg = pvt_region_idx_[well_cell_top];
if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) {
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells(), well_container.emplace_back(new StandardWell<TypeTag>(well_ecl,
param_, *rateConverter_, pvtreg, numComponents() ) ); time_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
w,
first_perf_index_[w],
well_perf_data_[w]));
} else { } else {
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells(), well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl,
param_, *rateConverter_, pvtreg, numComponents() ) ); time_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
w,
first_perf_index_[w],
well_perf_data_[w]));
} }
if (wellIsStopped) if (wellIsStopped)
well_container.back()->stopWell(); well_container.back()->stopWell();
@ -732,30 +744,32 @@ namespace Opm {
const Well& well_ecl = wells_ecl_[index_well_ecl]; 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 // 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]; const int pvtreg = pvt_region_idx_[well_cell_top];
if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) {
return WellInterfacePtr(new StandardWell<TypeTag>(well_ecl, report_step, wells(), return WellInterfacePtr(new StandardWell<TypeTag>(well_ecl,
param_, *rateConverter_, pvtreg, numComponents() ) ); report_step,
param_,
*rateConverter_,
pvtreg,
numComponents(),
numPhases(),
index_well_ecl,
first_perf_index_[index_well_ecl],
well_perf_data_[index_well_ecl]));
} else { } else {
return WellInterfacePtr(new MultisegmentWell<TypeTag>(well_ecl, report_step, wells(), return WellInterfacePtr(new MultisegmentWell<TypeTag>(well_ecl,
param_, *rateConverter_, pvtreg, numComponents() ) ); 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_) { for (const auto& well : well_container_) {
const bool needed_for_summary = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) || const bool needed_for_summary = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) ||
summaryConfig.hasSummaryKey( "WOPI:" + 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( "WWPP:" + well->name()) ||
summaryConfig.hasSummaryKey( "WOPP:" + 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(); const Well& eclWell = well->wellEcl();
bool needPotentialsForGuideRate = eclWell.getGuideRatePhase() == Well::GuideRateTarget::UNDEFINED; bool needPotentialsForGuideRate = eclWell.getGuideRatePhase() == Well::GuideRateTarget::UNDEFINED;
@ -1411,14 +1425,14 @@ namespace Opm {
int int
BlackoilWellModel<TypeTag>:: numLocalWells() const BlackoilWellModel<TypeTag>:: numLocalWells() const
{ {
return wells() ? wells()->number_of_wells : 0; return wells_ecl_.size();
} }
template<typename TypeTag> template<typename TypeTag>
int int
BlackoilWellModel<TypeTag>::numPhases() const BlackoilWellModel<TypeTag>::numPhases() const
{ {
return wells() ? wells()->number_of_phases : 1; return phase_usage_.num_phases;
} }
template<typename TypeTag> template<typename TypeTag>
@ -1560,22 +1574,14 @@ namespace Opm {
template<typename TypeTag> template<typename TypeTag>
bool bool
BlackoilWellModel<TypeTag>:: BlackoilWellModel<TypeTag>::
anyMSWellOpenLocal(const Wells* wells) const anyMSWellOpenLocal() const
{ {
bool any_ms_well_open = false; for (const auto& well : wells_ecl_) {
if (well.isMultiSegment()) {
const int nw = wells->number_of_wells; return true;
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;
} }
} }
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. // TODO: for now, we only use one type to save some implementation efforts, while improve later.
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell; 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 ModelParameters& param,
const RateConverterType& rate_converter, const RateConverterType& rate_converter,
const int pvtRegionIdx, 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, virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
@ -185,7 +189,6 @@ namespace Opm
using Base::well_cells_; using Base::well_cells_;
using Base::param_; using Base::param_;
using Base::well_index_; using Base::well_index_;
using Base::well_type_;
using Base::first_perf_; using Base::first_perf_;
using Base::saturation_table_number_; using Base::saturation_table_number_;
using Base::well_efficiency_factor_; 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 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 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; void assemblePressureEq(const int seg) const;
// hytrostatic pressure loss // hytrostatic pressure loss

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

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

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

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

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

@ -118,11 +118,15 @@ namespace Opm
compositionSwitchEnabled, compositionSwitchEnabled,
Indices::numPhases >; Indices::numPhases >;
/// Constructor /// Constructor
WellInterface(const Well& well, const int time_step, const Wells* wells, WellInterface(const Well& well, const int time_step,
const ModelParameters& param, const ModelParameters& param,
const RateConverterType& rate_converter, const RateConverterType& rate_converter,
const int pvtRegionIdx, 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 /// Virutal destructor
virtual ~WellInterface() {} virtual ~WellInterface() {}
@ -130,15 +134,18 @@ namespace Opm
/// Well name. /// Well name.
const std::string& name() const; 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 /// Index of well in the wells struct and wellState
int indexOfWell() const; int indexOfWell() const;
/// Well cells. /// Well cells.
const std::vector<int>& cells() const {return well_cells_; } const std::vector<int>& cells() const {return well_cells_; }
/// Well type, INJECTOR or PRODUCER.
WellType wellType() const;
/// Well controls /// Well controls
WellControls* wellControls() const; WellControls* wellControls() const;
@ -286,30 +293,12 @@ namespace Opm
const int current_step_; const int current_step_;
// the index of well in Wells struct
int index_of_well_;
// simulation parameters // simulation parameters
const ModelParameters& param_; 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 // number of the perforations for this well
int number_of_perforations_; int number_of_perforations_;
// record the index of the first perforation
// of states of individual well.
int first_perf_;
// well index for each perforation // well index for each perforation
std::vector<double> well_index_; std::vector<double> well_index_;
@ -372,6 +361,18 @@ namespace Opm
const int num_components_; 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_; std::vector<RateVector> connectionRates_;
bool wellIsStopped_; bool wellIsStopped_;

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

@ -27,75 +27,58 @@ namespace Opm
template<typename TypeTag> template<typename TypeTag>
WellInterface<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 ModelParameters& param,
const RateConverterType& rate_converter, const RateConverterType& rate_converter,
const int pvtRegionIdx, 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) : well_ecl_(well)
, current_step_(time_step) , current_step_(time_step)
, param_(param) , param_(param)
, rateConverter_(rate_converter) , rateConverter_(rate_converter)
, pvtRegionIdx_(pvtRegionIdx) , pvtRegionIdx_(pvtRegionIdx)
, num_components_(num_components) , num_components_(num_components)
, number_of_phases_(num_phases)
, index_of_well_(index_of_well)
, first_perf_(first_perf_index)
{ {
if (time_step < 0) { if (time_step < 0) {
OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface"); OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
} }
if (!wells) { ref_depth_ = well.getRefDepth();
OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface");
}
const std::string& well_name = well.name(); // We do not want to count SHUT perforations here, so
// it would be wrong to use wells.getConnections().size().
// looking for the location of the well in the wells struct number_of_perforations_ = perf_data.size();
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];
// perforations related // 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_); 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_); 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_); saturation_table_number_.resize(number_of_perforations_);
std::copy(wells->sat_table_id + perf_index_begin, int perf = 0;
wells->sat_table_id + perf_index_end, for (const auto& pd : perf_data) {
saturation_table_number_.begin() ); 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 // initialization of the completions mapping
@ -205,11 +188,22 @@ namespace Opm
template<typename TypeTag> template<typename TypeTag>
WellType bool
WellInterface<TypeTag>:: 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 // currently, we only updateWellTestState for producers
if (wellType() != PRODUCER) { if (this->isInjector()) {
return; return;
} }
@ -930,7 +924,8 @@ namespace Opm
bool allCompletionsClosed = true; bool allCompletionsClosed = true;
const auto& connections = well_ecl_.getConnections(); const auto& connections = well_ecl_.getConnections();
for (const auto& connection : connections) { 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; 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. // 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 // 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. // 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(); const int table_id = well_ecl_.vfp_table_number();
if (table_id <= 0) { if (table_id <= 0) {
return false; return false;
@ -1264,12 +1259,14 @@ namespace Opm
const auto& connections = well_ecl_.getConnections(); const auto& connections = well_ecl_.getConnections();
int perfIdx = 0; int perfIdx = 0;
for (const auto& connection : connections) { for (const auto& connection : connections) {
if (connection.state() == Connection::State::OPEN) {
if (wellTestState.hasCompletion(name(), connection.complnum())) { if (wellTestState.hasCompletion(name(), connection.complnum())) {
well_index_[perfIdx] = 0.0; well_index_[perfIdx] = 0.0;
} }
perfIdx++; perfIdx++;
} }
} }
}
template<typename TypeTag> template<typename TypeTag>
void void
@ -1295,7 +1292,7 @@ namespace Opm
void void
WellInterface<TypeTag>::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) 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 (well_ecl_.getDrainageRadius() < 0) {
if (new_well && perfIdx == 0) { if (new_well && perfIdx == 0) {
deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero"); deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero");

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

@ -24,6 +24,8 @@
#include <opm/core/wells.h> #include <opm/core/wells.h>
#include <opm/core/well_controls.h> #include <opm/core/well_controls.h>
#include <opm/output/data/Wells.hpp> #include <opm/output/data/Wells.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
#include <array> #include <array>
#include <map> #include <map>
@ -42,139 +44,54 @@ namespace Opm
typedef std::array< int, 3 > mapentry_t; typedef std::array< int, 3 > mapentry_t;
typedef std::map< std::string, mapentry_t > WellMapType; 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. /// Allocate and initialize if wells is non-null.
/// Also tries to give useful initial values to the bhp() and /// Also tries to give useful initial values to the bhp() and
/// wellRates() fields, depending on controls. The /// wellRates() fields, depending on controls. The
/// perfRates() field is filled with zero, and perfPress() /// perfRates() field is filled with zero, and perfPress()
/// with -1e100. /// 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 // clear old name mapping
wellMap_.clear(); wellMap_.clear();
wells_.reset( clone_wells( wells ) );
if (wells) { well_perf_data_ = well_perf_data;
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases; {
bhp_.resize(nw); // const int nw = wells->number_of_wells;
thp_.resize(nw); 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 temperature_.resize(nw, 273.15 + 20); // standard temperature for now
wellrates_.resize(nw * np, 0.0); wellrates_.resize(nw * np, 0.0);
int connpos = 0;
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER)); const Well& well = wells_ecl[w];
const WellControls* ctrl = wells->ctrls[w];
const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
// setup wellname -> well index mapping // Initialize bhp(), thp(), wellRates().
{ initSingleWell(cellPressures, w, well, pu, summary_state);
assert( wells->name[ w ] );
std::string name( wells->name[ w ] ); // 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 ); assert( name.size() > 0 );
mapentry_t& wellMapEntry = wellMap_[name]; mapentry_t& wellMapEntry = wellMap_[name];
wellMapEntry[ 0 ] = w; wellMapEntry[ 0 ] = w;
wellMapEntry[ 1 ] = wells->well_connpos[w]; wellMapEntry[ 1 ] = connpos;
// also store the number of perforations in this well
wellMapEntry[ 2 ] = num_perf_this_well; wellMapEntry[ 2 ] = num_perf_this_well;
} connpos += num_perf_this_well;
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;
}
}
} }
// The perforation rates and perforation pressures are // The perforation rates and perforation pressures are
// not expected to be consistent with bhp_ and wellrates_ // not expected to be consistent with bhp_ and wellrates_
// after init(). // after init().
perfrates_.resize(wells->well_connpos[nw], 0.0); perfrates_.resize(connpos, 0.0);
perfpress_.resize(wells->well_connpos[nw], -1e100); perfpress_.resize(connpos, -1e100);
} }
} }
@ -264,14 +181,11 @@ namespace Opm
well.rates.set( rt::gas, wv[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ] ); 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 ] const int num_perf_well = this->well_perf_data_[well_index].size();
- this->wells_->well_connpos[ well_index ];
well.connections.resize(num_perf_well); well.connections.resize(num_perf_well);
for( int i = 0; i < num_perf_well; ++i ) { for( int i = 0; i < num_perf_well; ++i ) {
const auto wi = this->wells_->well_connpos[ well_index ] + i; const auto active_index = this->well_perf_data_[well_index][i].cell_index;
const auto active_index = this->wells_->well_cells[ wi ];
auto& connection = well.connections[ i ]; auto& connection = well.connections[ i ];
connection.index = globalCellIdxMap[active_index]; connection.index = globalCellIdxMap[active_index];
connection.pressure = this->perfPress()[ itr.second[1] + i ]; connection.pressure = this->perfPress()[ itr.second[1] + i ];
@ -284,32 +198,10 @@ namespace Opm
} }
virtual ~WellState() {} virtual ~WellState() = default;
WellState() = default; WellState() = default;
WellState( const WellState& rhs ) : WellState(const WellState& rhs) = default;
bhp_( rhs.bhp_ ), WellState& operator=(const WellState& rhs) = default;
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;
}
private: private:
std::vector<double> bhp_; std::vector<double> bhp_;
@ -321,11 +213,139 @@ namespace Opm
WellMapType wellMap_; 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: protected:
struct wdel { std::vector<std::vector<PerforationData>> well_perf_data_;
void operator()( Wells* w ) { destroy_wells( w ); }
};
std::unique_ptr< Wells, wdel > wells_;
}; };
} // namespace Opm } // namespace Opm

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

@ -52,6 +52,8 @@ namespace Opm
public: public:
typedef BaseType :: WellMapType WellMapType; typedef BaseType :: WellMapType WellMapType;
virtual ~WellStateFullyImplicitBlackoil() = default;
// TODO: same definition with WellInterface, eventually they should go to a common header file. // TODO: same definition with WellInterface, eventually they should go to a common header file.
static const int Water = BlackoilPhases::Aqua; static const int Water = BlackoilPhases::Aqua;
static const int Oil = BlackoilPhases::Liquid; static const int Oil = BlackoilPhases::Liquid;
@ -67,29 +69,29 @@ namespace Opm
/// Allocate and initialize if wells is non-null. Also tries /// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates() /// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls /// and perfPhaseRates() fields, depending on controls
void init(const Wells* wells, void init(const std::vector<double>& cellPressures,
const std::vector<double>& cellPressures,
const Schedule& schedule, const Schedule& schedule,
const std::vector<Well>& wells_ecl, const std::vector<Well>& wells_ecl,
const int report_step, const int report_step,
const WellStateFullyImplicitBlackoil* prevState, 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 // 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 const int nw = wells_ecl.size();
if (wells == 0) {
return;
}
const int nw = wells->number_of_wells;
if( nw == 0 ) return ; if( nw == 0 ) return ;
// Initialize perfphaserates_, which must be done here. // Initialize perfphaserates_, which must be done here.
const int np = wells->number_of_phases; const int np = pu.num_phases;
const int nperf = wells->well_connpos[nw];
int nperf = 0;
for (const auto& wpd : well_perf_data) {
nperf += wpd.size();
}
well_reservoir_rates_.resize(nw * np, 0.0); well_reservoir_rates_.resize(nw * np, 0.0);
well_dissolved_gas_rates_.resize(nw, 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 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
+ ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::PRODUCTION_UPDATE
+ ScheduleEvents::INJECTION_UPDATE; + ScheduleEvents::INJECTION_UPDATE;
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
const int nw_wells_ecl = wells_ecl.size(); effective_events_occurred_[w]
int index_well_ecl = 0; = schedule.hasWellEvent(wells_ecl[w].name(), effective_events_mask, report_step);
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) );
} }
} // end of if (!well_ecl.empty() ) } // end of if (!well_ecl.empty() )
@ -139,21 +127,20 @@ namespace Opm
perf_skin_pressure_.clear(); perf_skin_pressure_.clear();
perf_skin_pressure_.resize(nperf, 0.0); perf_skin_pressure_.resize(nperf, 0.0);
int connpos = 0;
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER)); // Initialize perfphaserates_ to well
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
// rates divided by the number of perforations. // rates divided by the number of perforations.
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { const int num_perf_this_well = well_perf_data[w].size();
if (well_controls_well_is_open(ctrl)) { 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) { for (int p = 0; p < np; ++p) {
perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well); 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); current_injection_controls_.resize(nw);
@ -167,11 +154,12 @@ namespace Opm
// intialize wells that have been there before // intialize wells that have been there before
// order may change so the mapping is based on the well name // order may change so the mapping is based on the well name
if (prevState && !prevState->wellMap().empty()) { if (prevState && !prevState->wellMap().empty()) {
typedef typename WellMapType :: const_iterator const_iterator; connpos = 0;
const_iterator end = prevState->wellMap().end(); auto end = prevState->wellMap().end();
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
const std::string name( wells->name[ w ] ); const Well& well = wells_ecl[w];
const_iterator it = prevState->wellMap().find( name ); const int num_perf_this_well = well_perf_data[w].size();
auto it = prevState->wellMap().find(well.name());
if ( it != end ) if ( it != end )
{ {
const int oldIndex = (*it).second[ 0 ]; const int oldIndex = (*it).second[ 0 ];
@ -210,19 +198,18 @@ namespace Opm
// perfPhaseRates // perfPhaseRates
const int oldPerf_idx_beg = (*it).second[ 1 ]; const int oldPerf_idx_beg = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ]; 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, // copy perforation rates when the number of perforations is equal,
// otherwise initialize perfphaserates to well rates divided by the number of perforations. // otherwise initialize perfphaserates to well rates divided by the number of perforations.
if( num_perf_old_well == num_perf_this_well ) if( num_perf_old_well == num_perf_this_well )
{ {
int old_perf_phase_idx = oldPerf_idx_beg *np; int old_perf_phase_idx = oldPerf_idx_beg *np;
for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np; for (int perf_phase_idx = connpos*np;
perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx ) 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 ]; perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
} }
} else { } 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) { for (int p = 0; p < np; ++p) {
perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well); 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 ) if( num_perf_old_well == num_perf_this_well )
{ {
int oldPerf_idx = oldPerf_idx_beg; int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ]; for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{ {
perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ]; perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ];
} }
@ -243,8 +229,7 @@ namespace Opm
if( num_perf_old_well == num_perf_this_well ) if( num_perf_old_well == num_perf_this_well )
{ {
int oldPerf_idx = oldPerf_idx_beg; int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ]; for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{ {
perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ]; perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ];
} }
@ -258,8 +243,7 @@ namespace Opm
if( num_perf_old_well == num_perf_this_well ) if( num_perf_old_well == num_perf_this_well )
{ {
int oldPerf_idx = oldPerf_idx_beg; int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ]; for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx )
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{ {
perf_water_throughput_[ perf ] = prevState->perfThroughput()[ oldPerf_idx ]; perf_water_throughput_[ perf ] = prevState->perfThroughput()[ oldPerf_idx ];
perf_skin_pressure_[ perf ] = prevState->perfSkinPressure()[ 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 // If in the new step, there is no THP related target/limit anymore, its thp value should be
// set to zero. // set to zero.
const WellControls* ctrl = wells->ctrls[w]; const bool has_thp = well.isInjector() ? well.injectionControls(summary_state).hasControl(Well::InjectorCMode::THP)
const int nwc = well_controls_get_num(ctrl); : well.productionControls(summary_state).hasControl(Well::ProducerCMode::THP);
int ctrl_index = 0; if (!has_thp) {
for (; ctrl_index < nwc; ++ctrl_index) { thp()[w] = 0.0;
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
break;
}
}
// not finding any thp related control/limits
if (ctrl_index == nwc) {
thp()[w] = 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, void resize(const std::vector<Well>& wells_ecl,
const bool handle_ms_well, const size_t numCells, const PhaseUsage& pu) 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 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) { 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; 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 /* this is a reference or example on **how** to convert from
* WellState to something understood by opm-output. it is intended * WellState to something understood by opm-output. it is intended
* to be properly implemented and maintained as a part of * 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]] ); 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::dissolved_gas, this->well_dissolved_gas_rates_[w] );
well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_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) { for( auto& comp : well.connections) {
const auto rates = this->perfPhaseRates().begin() const auto rates = this->perfPhaseRates().begin()
+ (np * wt.second[ 1 ]) + (np * wt.second[ 1 ])
@ -505,7 +486,7 @@ namespace Opm
comp.rates.set( phs[ i ], *(rates + i) ); 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); const auto nseg = this->numSegments(w);
for (auto seg_ix = 0*nseg; seg_ix < nseg; ++seg_ix) { for (auto seg_ix = 0*nseg; seg_ix < nseg; ++seg_ix) {
@ -520,11 +501,11 @@ namespace Opm
/// init the MS well related. /// 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) const PhaseUsage& pu, const WellStateFullyImplicitBlackoil* prev_well_state)
{ {
// still using the order in wells // still using the order in wells
const int nw = wells->number_of_wells; const int nw = wells_ecl.size();
if (nw == 0) { if (nw == 0) {
return; return;
} }
@ -538,24 +519,12 @@ namespace Opm
seg_number_.clear(); seg_number_.clear();
nseg_ = 0; nseg_ = 0;
int connpos = 0;
// in the init function, the well rates and perforation rates have been initialized or copied from prevState // 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 // what we do here, is to set the segment rates and perforation rates
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
const int nw_wells_ecl = wells_ecl.size(); const auto& well_ecl = wells_ecl[w];
int index_well_ecl = 0; int num_perf_this_well = well_perf_data_[w].size();
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];
top_segment_index_.push_back(nseg_); top_segment_index_.push_back(nseg_);
if ( !well_ecl.isMultiSegment() ) { // not multi-segment well if ( !well_ecl.isMultiSegment() ) { // not multi-segment well
nseg_ += 1; nseg_ += 1;
@ -571,7 +540,6 @@ namespace Opm
const WellConnections& completion_set = well_ecl.getConnections(); const WellConnections& completion_set = well_ecl.getConnections();
// number of segment for this single well // number of segment for this single well
const int well_nseg = segment_set.size(); const int well_nseg = segment_set.size();
// const int nperf = completion_set.size();
int n_activeperf = 0; int n_activeperf = 0;
nseg_ += well_nseg; nseg_ += well_nseg;
for (auto segID = 0*well_nseg; segID < well_nseg; ++segID) { 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. // for the segrates_, now it becomes a recursive solution procedure.
{ {
const int np = numPhases(); const int np = numPhases();
const int start_perf = wells->well_connpos[w]; const int start_perf = connpos;
const int start_perf_next_well = wells->well_connpos[w + 1]; const int start_perf_next_well = connpos + num_perf_this_well;
// make sure the information from wells_ecl consistent with wells // make sure the information from wells_ecl consistent with wells
assert(n_activeperf == (start_perf_next_well - start_perf)); 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 // top segment is always the first one, and its pressure is the well bhp
segpress_.push_back(bhp()[w]); segpress_.push_back(bhp()[w]);
const int top_segment = top_segment_index_[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) { for (int seg = 1; seg < well_nseg; ++seg) {
if ( !segment_perforations[seg].empty() ) { if ( !segment_perforations[seg].empty() ) {
const int first_perf = segment_perforations[seg][0]; 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(segpress_.size()) == nseg_);
assert(int(segrates_.size()) == nseg_ * numPhases() ); assert(int(segrates_.size()) == nseg_ * numPhases() );
@ -663,8 +632,7 @@ namespace Opm
const auto& end = prev_well_state->wellMap().end(); const auto& end = prev_well_state->wellMap().end();
const int np = numPhases(); const int np = numPhases();
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
const std::string name( wells->name[w] ); const auto& it = prev_well_state->wellMap().find( wells_ecl[w].name() );
const auto& it = prev_well_state->wellMap().find( name );
if (it != end) { // the well is found in the prev_well_state 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 // 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 /// One rate pr well
double solventWellRate(const int w) const { 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; 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]; solvent_well_rate += perfRateSolvent_[perf];
} }
return solvent_well_rate; 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/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.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; using namespace Opm;
@ -42,54 +36,24 @@ BOOST_AUTO_TEST_CASE(TestStoppedWells)
Opm::Parser parser; Opm::Parser parser;
Opm::Deck deck(parser.parseFile(filename)); Opm::Deck deck(parser.parseFile(filename));
Opm::EclipseState eclipseState(deck); Opm::EclipseState eclipseState(deck);
Opm::GridManager vanguard(eclipseState.getInputGrid());
const auto& grid = eclipseState.getInputGrid(); const auto& grid = eclipseState.getInputGrid();
const TableManager table ( deck ); const TableManager table ( deck );
const Eclipse3DProperties eclipseProperties ( deck , table, grid); const Eclipse3DProperties eclipseProperties ( deck , table, grid);
const Opm::Runspec runspec (deck); const Opm::Runspec runspec (deck);
const Schedule sched(deck, grid, eclipseProperties, runspec); 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 // Both wells are open in the first schedule step
{ {
Opm::WellsManager wellsManager(eclipseState , sched, summaryState, 0, *vanguard.c_grid()); auto wells = sched.getWells(0);
const Wells* wells = wellsManager.c_wells(); BOOST_CHECK(wells[0].getStatus() == Opm::Well::Status::OPEN);
const struct WellControls* ctrls0 = wells->ctrls[0]; BOOST_CHECK(wells[1].getStatus() == Opm::Well::Status::OPEN);
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
} }
// The injector is stopped // The injector is stopped
{ {
Opm::WellsManager wellsManager(eclipseState, sched, summaryState, 1 , *vanguard.c_grid()); auto wells = sched.getWells(1);
const Wells* wells = wellsManager.c_wells(); BOOST_CHECK(wells[0].getStatus() == Opm::Well::Status::STOP);
const struct WellControls* ctrls0 = wells->ctrls[0]; BOOST_CHECK(wells[1].getStatus() == Opm::Well::Status::OPEN);
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
} }
} }

58
tests/test_wellmodel.cpp
View File

@ -77,36 +77,13 @@ struct SetupTest {
schedule.reset( new Opm::Schedule(deck, ecl_state->getInputGrid(), eclipseProperties, runspec)); 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()))); 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; 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::EclipseState> ecl_state;
std::unique_ptr<const Opm::Schedule> schedule; std::unique_ptr<const Opm::Schedule> schedule;
std::unique_ptr<Opm::SummaryState> summaryState; std::unique_ptr<Opm::SummaryState> summaryState;
std::vector<std::vector<Opm::PerforationData>> well_perf_data;
int current_timestep; int current_timestep;
}; };
@ -132,7 +109,6 @@ BOOST_GLOBAL_FIXTURE(GlobalFixture);
BOOST_AUTO_TEST_CASE(TestStandardWellInput) { BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const SetupTest setup_test; 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); const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep);
BOOST_CHECK_EQUAL( wells_ecl.size(), 2); BOOST_CHECK_EQUAL( wells_ecl.size(), 2);
const Opm::Well& well = wells_ecl[1]; const Opm::Well& well = wells_ecl[1];
@ -149,36 +125,24 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
rateConverter.reset(new RateConverterType (phaseUsage, rateConverter.reset(new RateConverterType (phaseUsage,
std::vector<int>(10, 0))); std::vector<int>(10, 0)));
const int pvtIdx = 0; Opm::PerforationData dummy;
const int num_comp = wells->number_of_phases; 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, -1, param, *rateConverter, 0, 3, 3, 0, 0, pdata), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr , param, *rateConverter, pvtIdx, num_comp), std::invalid_argument);
} }
BOOST_AUTO_TEST_CASE(TestBehavoir) { BOOST_AUTO_TEST_CASE(TestBehavoir) {
const SetupTest setup_test; 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 auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep);
const int current_timestep = setup_test.current_timestep; const int current_timestep = setup_test.current_timestep;
std::vector<std::unique_ptr<const StandardWell> > wells; 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; const Opm::BlackoilModelParametersEbos<TTAG(EclFlowProblem)> param;
for (int w = 0; w < nw; ++w) { 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 // For the conversion between the surface volume rate and resrevoir voidage rate
typedef Opm::BlackOilFluidSystem<double> FluidSystem; typedef Opm::BlackOilFluidSystem<double> FluidSystem;
using RateConverterType = Opm::RateConverter:: using RateConverterType = Opm::RateConverter::
@ -191,11 +155,9 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
// Compute reservoir volumes for RESV controls. // Compute reservoir volumes for RESV controls.
rateConverter.reset(new RateConverterType (phaseUsage, rateConverter.reset(new RateConverterType (phaseUsage,
std::vector<int>(10, 0))); std::vector<int>(10, 0)));
Opm::PerforationData dummy;
const int pvtIdx = 0; std::vector<Opm::PerforationData> pdata(wells_ecl[w].getConnections().size(), dummy);
const int num_comp = wells_struct->number_of_phases; wells.emplace_back(new StandardWell(wells_ecl[w], current_timestep, param, *rateConverter, 0, 3, 3, w, 0, pdata) );
wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param, *rateConverter, pvtIdx, num_comp) );
} }
} }
@ -203,7 +165,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{ {
const auto& well = wells[0]; const auto& well = wells[0];
BOOST_CHECK_EQUAL(well->name(), "PROD1"); BOOST_CHECK_EQUAL(well->name(), "PROD1");
BOOST_CHECK(well->wellType() == PRODUCER); BOOST_CHECK(well->isProducer());
BOOST_CHECK(well->numEq == 3); BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4); BOOST_CHECK(well->numStaticWellEq== 4);
} }
@ -212,7 +174,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{ {
const auto& well = wells[1]; const auto& well = wells[1];
BOOST_CHECK_EQUAL(well->name(), "INJE1"); BOOST_CHECK_EQUAL(well->name(), "INJE1");
BOOST_CHECK(well->wellType() == INJECTOR); BOOST_CHECK(well->isInjector());
BOOST_CHECK(well->numEq == 3); BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4); BOOST_CHECK(well->numStaticWellEq== 4);
} }

58
tests/test_wellstatefullyimplicitblackoil.cpp
View File

@ -53,13 +53,60 @@ struct Setup
, grid (es.getInputGrid()) , grid (es.getInputGrid())
, sched(deck, es) , sched(deck, es)
, st(std::chrono::system_clock::from_time_t(sched.getStartTime())) , 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::EclipseState es;
Opm::PhaseUsage pu; Opm::PhaseUsage pu;
Opm::GridManager grid; Opm::GridManager grid;
Opm::Schedule sched; Opm::Schedule sched;
Opm::SummaryState st; Opm::SummaryState st;
std::vector<std::vector<Opm::PerforationData>> well_perf_data;
}; };
namespace { namespace {
@ -72,14 +119,11 @@ namespace {
std::vector<double>(setup.grid.c_grid()->number_of_cells, std::vector<double>(setup.grid.c_grid()->number_of_cells,
100.0*Opm::unit::barsa); 100.0*Opm::unit::barsa);
const Opm::WellsManager wmgr{setup.es, setup.sched, setup.st, timeStep, *setup.grid.c_grid()}; state.init(cpress, setup.sched,
state.init(wmgr.c_wells(), cpress, setup.sched,
setup.sched.getWells(timeStep), setup.sched.getWells(timeStep),
timeStep, nullptr, setup.pu); timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st);
state.initWellStateMSWell(wmgr.c_wells(), state.initWellStateMSWell(setup.sched.getWells(timeStep),
setup.sched.getWells(timeStep),
setup.pu, nullptr); setup.pu, nullptr);
return state; return state;