From 87188f5862584c737c7d61e8ea7b29906df56b53 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 23 Oct 2019 09:09:45 +0200 Subject: [PATCH] Avoid using the Wells struct. --- CMakeLists_files.cmake | 1 + opm/simulators/wells/BlackoilWellModel.hpp | 13 +- .../wells/BlackoilWellModel_impl.hpp | 320 +++++++++--------- opm/simulators/wells/MultisegmentWell.hpp | 10 +- .../wells/MultisegmentWell_impl.hpp | 60 ++-- opm/simulators/wells/PerforationData.hpp | 36 ++ opm/simulators/wells/StandardWell.hpp | 11 +- opm/simulators/wells/StandardWell_impl.hpp | 216 +++++++----- opm/simulators/wells/WellInterface.hpp | 47 +-- opm/simulators/wells/WellInterface_impl.hpp | 113 +++---- opm/simulators/wells/WellState.hpp | 316 +++++++++-------- .../wells/WellStateFullyImplicitBlackoil.hpp | 177 ++++------ tests/test_stoppedwells.cpp | 48 +-- tests/test_wellmodel.cpp | 58 +--- tests/test_wellstatefullyimplicitblackoil.cpp | 58 +++- 15 files changed, 782 insertions(+), 702 deletions(-) create mode 100644 opm/simulators/wells/PerforationData.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index eb8761605..3f8c9904b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -187,6 +187,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/utils/gatherDeferredLogger.hpp opm/simulators/utils/moduleVersion.hpp opm/simulators/utils/ParallelRestart.hpp + opm/simulators/wells/PerforationData.hpp opm/simulators/wells/RateConverter.hpp opm/simulators/wells/SimFIBODetails.hpp opm/simulators/wells/WellConnectionAuxiliaryModule.hpp diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index b07e11651..77dde08a9 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -40,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -233,8 +234,10 @@ namespace Opm { protected: Simulator& ebosSimulator_; - std::unique_ptr wells_manager_; + std::vector< Well > wells_ecl_; + std::vector< std::vector > well_perf_data_; + std::vector first_perf_index_; bool wells_active_; @@ -246,8 +249,10 @@ namespace Opm { std::vector is_cell_perforated_; + void initializeWellPerfData(); + // create the well container - std::vector createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger); + std::vector createWellContainer(const int time_step); WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, Opm::DeferredLogger& deferred_logger) const; @@ -278,8 +283,6 @@ namespace Opm { // used to better efficiency of calcuation mutable BVector scaleAddRes_; - const Wells* wells() const { return wells_manager_->c_wells(); } - const Grid& grid() const { return ebosSimulator_.vanguard().grid(); } @@ -373,7 +376,7 @@ namespace Opm { WellStateFullyImplicitBlackoil& state ) const; // whether there exists any multisegment well open on this process - bool anyMSWellOpenLocal(const Wells* wells) const; + bool anyMSWellOpenLocal() const; const Well& getWellEcl(const std::string& well_name) const; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 861790258..f5ee5340b 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -37,6 +37,19 @@ namespace Opm { // Create the guide rate container. guideRate_.reset(new GuideRate (ebosSimulator_.vanguard().schedule())); + + // calculate the number of elements of the compressed sequential grid. this needs + // to be done in two steps because the dune communicator expects a reference as + // argument for sum() + const auto& gridView = ebosSimulator_.gridView(); + number_of_cells_ = gridView.size(/*codim=*/0); + global_nc_ = gridView.comm().sum(number_of_cells_); + + // Set up cartesian mapping. + const auto& grid = ebosSimulator_.vanguard().grid(); + const auto& cartDims = Opm::UgGridHelpers::cartDims(grid); + setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid), + cartDims[0]*cartDims[1]*cartDims[2]); } template @@ -46,31 +59,15 @@ namespace Opm { { const Opm::EclipseState& eclState = ebosSimulator_.vanguard().eclState(); - gravity_ = ebosSimulator_.problem().gravity()[2]; - extractLegacyCellPvtRegionIndex_(); extractLegacyDepth_(); phase_usage_ = phaseUsageFromDeck(eclState); - const auto& gridView = ebosSimulator_.gridView(); - - // calculate the number of elements of the compressed sequential grid. this needs - // to be done in two steps because the dune communicator expects a reference as - // argument for sum() - number_of_cells_ = gridView.size(/*codim=*/0); - global_nc_ = gridView.comm().sum(number_of_cells_); gravity_ = ebosSimulator_.problem().gravity()[2]; - extractLegacyCellPvtRegionIndex_(); - extractLegacyDepth_(); initial_step_ = true; - const auto& grid = ebosSimulator_.vanguard().grid(); - const auto& cartDims = Opm::UgGridHelpers::cartDims(grid); - setupCartesianToCompressed_(Opm::UgGridHelpers::globalCell(grid), - cartDims[0]*cartDims[1]*cartDims[2]); - // add the eWoms auxiliary module for the wells to the list ebosSimulator_.model().addAuxiliaryModule(this); @@ -210,24 +207,18 @@ namespace Opm { Opm::DeferredLogger local_deferredLogger; const Grid& grid = ebosSimulator_.vanguard().grid(); - const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames(); - const auto& eclState = ebosSimulator_.vanguard().eclState(); const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - wells_ecl_ = schedule().getWells(timeStepIdx); - - // Create wells and well state. - wells_manager_.reset( new WellsManager (eclState, - schedule(), - summaryState, - timeStepIdx, - Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::globalCell(grid), - Opm::UgGridHelpers::cartDims(grid), - Opm::UgGridHelpers::dimensions(grid), - Opm::UgGridHelpers::cell2Faces(grid), - Opm::UgGridHelpers::beginFaceCentroids(grid), - grid.comm().size() > 1, - defunct_well_names) ); + // Make wells_ecl_ contain only this partition's non-shut wells. + { + const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames(); + auto is_shut_or_defunct = [&defunct_well_names](const Well& well) { + return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end()); + }; + auto w = schedule().getWells(timeStepIdx); + w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end()); + wells_ecl_.swap(w); + } + initializeWellPerfData(); // Wells are active if they are active wells on at least // one process. @@ -269,35 +260,22 @@ namespace Opm { } cellPressures[cellIdx] = perf_pressure; } - well_state_.init(wells(), cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_); + well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_, well_perf_data_, summaryState); // handling MS well related - if (param_.use_multisegment_well_&& anyMSWellOpenLocal(wells())) { // if we use MultisegmentWell model - well_state_.initWellStateMSWell(wells(), wells_ecl_, phase_usage_, &previous_well_state_); + if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model + well_state_.initWellStateMSWell(wells_ecl_, phase_usage_, &previous_well_state_); } - const int nw = wells()->number_of_wells; + const int nw = wells_ecl_.size(); for (int w = 0; w name[w]); - for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) { - if (well_name == wells_ecl_[index_well_ecl].name()) { - break; - } - } - - // It should be able to find in wells_ecl. - if (index_well_ecl == nw_wells_ecl) { - OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl "); - } - const auto well = wells_ecl_[index_well_ecl]; + const auto& well = wells_ecl_[w]; const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE + ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::INJECTION_UPDATE + ScheduleEvents::NEW_WELL; - if(!schedule().hasWellEvent(well_name, effective_events_mask, timeStepIdx)) + if(!schedule().hasWellEvent(well.name(), effective_events_mask, timeStepIdx)) continue; if (well.isProducer()) { @@ -347,7 +325,7 @@ namespace Opm { wellTesting(reportStepIdx, simulationTime, local_deferredLogger); // create the well container - well_container_ = createWellContainer(reportStepIdx, local_deferredLogger); + well_container_ = createWellContainer(reportStepIdx); // do the initialization for all the wells // TODO: to see whether we can postpone of the intialization of the well containers to @@ -400,9 +378,16 @@ namespace Opm { //compute well guideRates const int np = numPhases(); for (const auto& well : well_container_) { - const double& oilpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Liquid]]; - const double& gaspot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Vapour]]; - const double& waterpot = well_state_.wellPotentials()[well->indexOfWell() * np + phase_usage_.phase_pos[BlackoilPhases::Aqua]]; + const auto wpot = well_state_.wellPotentials().data() + well->indexOfWell()*np; + const double oilpot = (phase_usage_.phase_used[BlackoilPhases::Liquid] > 0) + ? wpot[phase_usage_.phase_pos[BlackoilPhases::Liquid]] + : 0.0; + const double gaspot = (phase_usage_.phase_used[BlackoilPhases::Vapour] > 0) + ? wpot[phase_usage_.phase_pos[BlackoilPhases::Vapour]] + : 0.0; + const double waterpot = (phase_usage_.phase_used[BlackoilPhases::Aqua] > 0) + ? wpot[phase_usage_.phase_pos[BlackoilPhases::Aqua]] + : 0.0; guideRate_->compute(well->name(), reportStepIdx, simulationTime, oilpot, gaspot, waterpot); } const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); @@ -470,7 +455,7 @@ namespace Opm { Opm::DeferredLogger local_deferredLogger; for (const auto& well : well_container_) { - if (GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well->wellType() == INJECTOR) { + if (GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well->isInjector()) { well->updateWaterThroughput(dt, well_state_); } } @@ -538,8 +523,6 @@ namespace Opm { BlackoilWellModel:: initFromRestartFile(const RestartValue& restartValues) { - const auto& defunctWellNames = ebosSimulator_.vanguard().defunctWellNames(); - // The restart step value is used to identify wells present at the given // time step. Wells that are added at the same time step as RESTART is initiated // will not be present in a restart file. Use the previous time step to retrieve @@ -547,29 +530,25 @@ namespace Opm { const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0); const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - WellsManager wellsmanager(eclState(), - schedule(), - summaryState, - report_step, - Opm::UgGridHelpers::numCells(grid()), - Opm::UgGridHelpers::globalCell(grid()), - Opm::UgGridHelpers::cartDims(grid()), - Opm::UgGridHelpers::dimensions(grid()), - Opm::UgGridHelpers::cell2Faces(grid()), - Opm::UgGridHelpers::beginFaceCentroids(grid()), - grid().comm().size() > 1, - defunctWellNames); + // Make wells_ecl_ contain only this partition's non-shut wells. + { + const auto& defunct_well_names = ebosSimulator_.vanguard().defunctWellNames(); + auto is_shut_or_defunct = [&defunct_well_names](const Well& well) { + return (well.getStatus() == Well::Status::SHUT) || (defunct_well_names.find(well.name()) != defunct_well_names.end()); + }; + auto w = schedule().getWells(report_step); + w.erase(std::remove_if(w.begin(), w.end(), is_shut_or_defunct), w.end()); + wells_ecl_.swap(w); + } - const Wells* wells = wellsmanager.c_wells(); + initializeWellPerfData(); - wells_ecl_ = schedule().getWells(report_step); - - const int nw = wells->number_of_wells; + const int nw = wells_ecl_.size(); if (nw > 0) { const auto phaseUsage = phaseUsageFromDeck(eclState()); const size_t numCells = Opm::UgGridHelpers::numCells(grid()); - const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal(wells)); - well_state_.resize(wells, wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage); // Resize for restart step + const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal()); + well_state_.resize(wells_ecl_, schedule(), handle_ms_well, numCells, phaseUsage, well_perf_data_, summaryState); // Resize for restart step wellsToState(restartValues.wells, phaseUsage, handle_ms_well, well_state_); } @@ -578,20 +557,7 @@ namespace Opm { const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST(); if (true || ecl_compatible_rst) { // always set the control from the schedule for (int w = 0; w name[w]); - for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) { - if (well_name == wells_ecl_[index_well_ecl].name()) { - break; - } - } - - // It should be able to find in wells_ecl. - if (index_well_ecl == nw_wells_ecl) { - OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl "); - } - const auto& well = wells_ecl_[index_well_ecl]; + const auto& well = wells_ecl_[w]; if (well.isProducer()) { const auto controls = well.productionControls(summaryState); @@ -613,10 +579,59 @@ namespace Opm { + template + void + BlackoilWellModel:: + 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 std::vector::WellInterfacePtr > BlackoilWellModel:: - createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger) + createWellContainer(const int time_step) { std::vector well_container; @@ -624,27 +639,9 @@ namespace Opm { if (nw > 0) { well_container.reserve(nw); - - // With the following way, it will have the same order with wells struct - // Hopefully, it can generate the same residual history with master branch for (int w = 0; w < nw; ++w) { - const std::string well_name = std::string(wells()->name[w]); - - // finding the location of the well in wells_ecl - const int nw_wells_ecl = wells_ecl_.size(); - int index_well = 0; - for (; index_well < nw_wells_ecl; ++index_well) { - if (well_name == wells_ecl_[index_well].name()) { - break; - } - } - - // It should be able to find in wells_ecl. - if (index_well == nw_wells_ecl) { - OPM_DEFLOG_THROW(std::runtime_error, "Could not find well " + well_name + " in wells_ecl ", deferred_logger); - } - - const Well& well_ecl = wells_ecl_[index_well]; + const Well& well_ecl = wells_ecl_[w]; + const std::string& well_name = well_ecl.name(); // A new WCON keywords can re-open a well that was closed/shut due to Physical limit if ( wellTestState_.hasWellClosed(well_name)) { @@ -663,7 +660,6 @@ namespace Opm { } } - // TODO: should we do this for all kinds of closing reasons? // something like wellTestState_.hasWell(well_name)? bool wellIsStopped = false; @@ -688,15 +684,31 @@ namespace Opm { } // Use the pvtRegionIdx from the top cell - const int well_cell_top = wells()->well_cells[wells()->well_connpos[w]]; + const int well_cell_top = well_perf_data_[w][0].cell_index; const int pvtreg = pvt_region_idx_[well_cell_top]; - if ( !well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { - well_container.emplace_back(new StandardWell(well_ecl, time_step, wells(), - param_, *rateConverter_, pvtreg, numComponents() ) ); + if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { + well_container.emplace_back(new StandardWell(well_ecl, + time_step, + param_, + *rateConverter_, + pvtreg, + numComponents(), + numPhases(), + w, + first_perf_index_[w], + well_perf_data_[w])); } else { - well_container.emplace_back(new MultisegmentWell(well_ecl, time_step, wells(), - param_, *rateConverter_, pvtreg, numComponents() ) ); + well_container.emplace_back(new MultisegmentWell(well_ecl, + time_step, + param_, + *rateConverter_, + pvtreg, + numComponents(), + numPhases(), + w, + first_perf_index_[w], + well_perf_data_[w])); } if (wellIsStopped) well_container.back()->stopWell(); @@ -732,30 +744,32 @@ namespace Opm { const Well& well_ecl = wells_ecl_[index_well_ecl]; - // Finding the location of the well in wells struct. - const int nw = numLocalWells(); - int well_index_wells = -999; - for (int w = 0; w < nw; ++w) { - if (well_name == std::string(wells()->name[w])) { - well_index_wells = w; - break; - } - } - - if (well_index_wells < 0) { - OPM_DEFLOG_THROW(std::logic_error, "Could not find the well " << well_name << " in the well struct ", deferred_logger); - } - // Use the pvtRegionIdx from the top cell - const int well_cell_top = wells()->well_cells[wells()->well_connpos[well_index_wells]]; + const int well_cell_top = well_perf_data_[index_well_ecl][0].cell_index; const int pvtreg = pvt_region_idx_[well_cell_top]; - if ( !well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { - return WellInterfacePtr(new StandardWell(well_ecl, report_step, wells(), - param_, *rateConverter_, pvtreg, numComponents() ) ); + if (!well_ecl.isMultiSegment() || !param_.use_multisegment_well_) { + return WellInterfacePtr(new StandardWell(well_ecl, + report_step, + param_, + *rateConverter_, + pvtreg, + numComponents(), + numPhases(), + index_well_ecl, + first_perf_index_[index_well_ecl], + well_perf_data_[index_well_ecl])); } else { - return WellInterfacePtr(new MultisegmentWell(well_ecl, report_step, wells(), - param_, *rateConverter_, pvtreg, numComponents() ) ); + return WellInterfacePtr(new MultisegmentWell(well_ecl, + report_step, + param_, + *rateConverter_, + pvtreg, + numComponents(), + numPhases(), + index_well_ecl, + first_perf_index_[index_well_ecl], + well_perf_data_[index_well_ecl])); } } @@ -1183,10 +1197,10 @@ namespace Opm { for (const auto& well : well_container_) { const bool needed_for_summary = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) || summaryConfig.hasSummaryKey( "WOPI:" + well->name()) || - summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->wellType() == INJECTOR) || + summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->isInjector()) || ((summaryConfig.hasSummaryKey( "WWPP:" + well->name()) || summaryConfig.hasSummaryKey( "WOPP:" + well->name()) || - summaryConfig.hasSummaryKey( "WGPP:" + well->name())) && well->wellType() == PRODUCER); + summaryConfig.hasSummaryKey( "WGPP:" + well->name())) && well->isProducer()); const Well& eclWell = well->wellEcl(); bool needPotentialsForGuideRate = eclWell.getGuideRatePhase() == Well::GuideRateTarget::UNDEFINED; @@ -1411,14 +1425,14 @@ namespace Opm { int BlackoilWellModel:: numLocalWells() const { - return wells() ? wells()->number_of_wells : 0; + return wells_ecl_.size(); } template int - BlackoilWellModel:: numPhases() const + BlackoilWellModel::numPhases() const { - return wells() ? wells()->number_of_phases : 1; + return phase_usage_.num_phases; } template @@ -1560,22 +1574,14 @@ namespace Opm { template bool BlackoilWellModel:: - anyMSWellOpenLocal(const Wells* wells) const + anyMSWellOpenLocal() const { - bool any_ms_well_open = false; - - const int nw = wells->number_of_wells; - for (int w = 0; w < nw; ++w) { - const std::string well_name = std::string(wells->name[w]); - - const Well& well_ecl = getWellEcl(well_name); - - if (well_ecl.isMultiSegment() ) { - any_ms_well_open = true; - break; + for (const auto& well : wells_ecl_) { + if (well.isMultiSegment()) { + return true; } } - return any_ms_well_open; + return false; } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 9af6f47b3..ae46ae645 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -98,11 +98,15 @@ namespace Opm // TODO: for now, we only use one type to save some implementation efforts, while improve later. typedef DenseAd::Evaluation EvalWell; - MultisegmentWell(const Well& well, const int time_step, const Wells* wells, + MultisegmentWell(const Well& well, const int time_step, const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, - const int num_components); + const int num_components, + const int num_phases, + const int index_of_well, + const int first_perf_index, + const std::vector& perf_data); virtual void init(const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, @@ -185,7 +189,6 @@ namespace Opm using Base::well_cells_; using Base::param_; using Base::well_index_; - using Base::well_type_; using Base::first_perf_; using Base::saturation_table_number_; using Base::well_efficiency_factor_; @@ -355,7 +358,6 @@ namespace Opm void assembleGroupProductionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger); void assembleGroupInjectionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, const Well::InjectorType& injectorType, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger); - void assemblePressureEq(const int seg) const; // hytrostatic pressure loss diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 4587ff827..27360f9fe 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -28,12 +28,16 @@ namespace Opm template MultisegmentWell:: - MultisegmentWell(const Well& well, const int time_step, const Wells* wells, + MultisegmentWell(const Well& well, const int time_step, const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, - const int num_components) - : Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components) + const int num_components, + const int num_phases, + const int index_of_well, + const int first_perf_index, + const std::vector& perf_data) + : Base(well, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, first_perf_index, perf_data) , segment_perforations_(numberOfSegments()) , segment_inlets_(numberOfSegments()) , cell_perforation_depth_diffs_(number_of_perforations_, 0.0) @@ -756,7 +760,7 @@ namespace Opm primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate; } } else { // total_seg_rate == 0 - if (well_type_ == INJECTOR) { + if (this->isInjector()) { // only single phase injection handled auto phase = well.getInjectionProperties().injectorType; @@ -776,7 +780,7 @@ namespace Opm } } - } else if (well_type_ == PRODUCER) { // producers + } else if (this->isProducer()) { // producers if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { primary_variables_[seg][WFrac] = 1.0 / number_of_phases_; } @@ -942,7 +946,7 @@ namespace Opm // make sure that no injector produce and no producer inject if (seg == 0) { - if (well_type_ == INJECTOR) { + if (this->isInjector()) { primary_variables_[seg][GTotal] = std::max( primary_variables_[seg][GTotal], 0.0); } else { primary_variables_[seg][GTotal] = std::min( primary_variables_[seg][GTotal], 0.0); @@ -1186,7 +1190,7 @@ namespace Opm // producing perforations if ( drawdown > 0.0) { // Do nothing is crossflow is not allowed - if (!allow_cf && well_type_ == INJECTOR) { + if (!allow_cf && this->isInjector()) { return; } @@ -1206,7 +1210,7 @@ namespace Opm } } else { // injecting perforations // Do nothing if crossflow is not allowed - if (!allow_cf && well_type_ == PRODUCER) { + if (!allow_cf && this->isProducer()) { return; } @@ -1263,7 +1267,7 @@ namespace Opm } // end for injection perforations // calculating the perforation solution gas rate and solution oil rates - if (well_type_ == PRODUCER) { + if (this->isProducer()) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); @@ -1516,7 +1520,7 @@ namespace Opm // the result will contain the derivative with resepct to GTotal in segment seg, // and the derivatives with respect to WFrac GFrac in segment seg_upwind. // the derivative with respect to SPres should be zero. - if (seg == 0 && well_type_ == INJECTOR) { + if (seg == 0 && this->isInjector()) { const Well& well = Base::wellEcl(); auto phase = well.getInjectionProperties().injectorType; @@ -1630,7 +1634,7 @@ namespace Opm if (wellIsStopped_) { control_eq = getSegmentGTotal(0); - } else if (well.isInjector() ) { + } else if (this->isInjector() ) { const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index]; const auto controls = well.injectionControls(summaryState); @@ -1711,13 +1715,16 @@ namespace Opm rates[ Gas ] = getSegmentRate(0, Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)); } control_eq = getSegmentPressure(0) - calculateBhpFromThp(rates, well, summaryState, deferred_logger); - break; } + break; + } + case Well::InjectorCMode::BHP: { const auto& bhp = controls.bhp_limit; control_eq = getSegmentPressure(0) - bhp; break; } + case Well::InjectorCMode::GRUP: { assert(well.isAvailableForGroupControl()); @@ -1725,6 +1732,7 @@ namespace Opm assembleGroupInjectionControl(group, well_state, schedule, summaryState, controls.injector_type, control_eq, efficiencyFactor, deferred_logger); break; } + case Well::InjectorCMode::CMODE_UNDEFINED: { OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + name(), deferred_logger); @@ -1823,7 +1831,8 @@ namespace Opm rates[ Gas ] = getSegmentRate(0, Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)); } control_eq = getSegmentPressure(0) - calculateBhpFromThp(rates, well, summaryState, deferred_logger); - break; } + break; + } case Well::ProducerCMode::GRUP: { assert(well.isAvailableForGroupControl()); @@ -1844,7 +1853,6 @@ namespace Opm } } - // using control_eq to update the matrix and residuals resWell_[0][SPres] = control_eq.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { @@ -1903,14 +1911,14 @@ namespace Opm const double rho = segment_densities_[0].value(); double thp = 0.0; - if (well_type_ == INJECTOR) { + if (this->isInjector()) { const int table_id = well_ecl_.vfp_table_number(); const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); } - else if (well_type_ == PRODUCER) { + else if (this->isProducer()) { const int table_id = well_ecl_.vfp_table_number(); const double alq = well_ecl_.alq_value(); const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); @@ -2675,7 +2683,7 @@ namespace Opm computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); // updating the solution gas rate and solution oil rate - if (well_type_ == PRODUCER) { + if (this->isProducer()) { well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate; well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; } @@ -2762,8 +2770,8 @@ namespace Opm // for now, if there is one perforation can produce/inject in the correct // direction, we consider this well can still produce/inject. // TODO: it can be more complicated than this to cause wrong-signed rates - if ( (drawdown < 0. && well_type_ == INJECTOR) || - (drawdown > 0. && well_type_ == PRODUCER) ) { + if ( (drawdown < 0. && this->isInjector()) || + (drawdown > 0. && this->isProducer()) ) { all_drawdown_wrong_direction = false; break; } @@ -3056,9 +3064,8 @@ namespace Opm { double control_tolerance = 0.; - const auto& well = well_ecl_; const int well_index = index_of_well_; - if (well.isInjector() ) + if (this->isInjector() ) { const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index]; switch(current) { @@ -3080,7 +3087,7 @@ namespace Opm } } - if (well.isProducer() ) + if (this->isProducer() ) { const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index]; switch(current) { @@ -3124,9 +3131,8 @@ namespace Opm using CR = ConvergenceReport; CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid; - const auto& well = well_ecl_; const int well_index = index_of_well_; - if (well.isInjector() ) + if (this->isInjector() ) { const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index]; switch(current) { @@ -3152,7 +3158,7 @@ namespace Opm } } - if (well.isProducer() ) + if (this->isProducer() ) { const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index]; switch(current) { @@ -3208,8 +3214,8 @@ namespace Opm // special treatment is needed for segment 0 if (seg == 0) { // we are not supposed to have injecting producers and producing injectors - assert( ! (well_type_ == PRODUCER && primary_variables_evaluation_[seg][GTotal] > 0.) ); - assert( ! (well_type_ == INJECTOR && primary_variables_evaluation_[seg][GTotal] < 0.) ); + assert( ! (this->isProducer() && primary_variables_evaluation_[seg][GTotal] > 0.) ); + assert( ! (this->isInjector() && primary_variables_evaluation_[seg][GTotal] < 0.) ); upwinding_segments_[seg] = seg; continue; } diff --git a/opm/simulators/wells/PerforationData.hpp b/opm/simulators/wells/PerforationData.hpp new file mode 100644 index 000000000..47438026d --- /dev/null +++ b/opm/simulators/wells/PerforationData.hpp @@ -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 . +*/ + +#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 diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 310c88be5..6022d593d 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -140,11 +140,15 @@ namespace Opm using Base::contiFoamEqIdx; static const int contiEnergyEqIdx = Indices::contiEnergyEqIdx; - StandardWell(const Well& well, const int time_step, const Wells* wells, + StandardWell(const Well& well, const int time_step, const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, - const int num_components); + const int num_components, + const int num_phases, + const int index_of_well, + const int first_perf_index, + const std::vector& perf_data); virtual void init(const PhaseUsage* phase_usage_arg, const std::vector& depth_arg, @@ -231,10 +235,8 @@ namespace Opm using Base::number_of_perforations_; using Base::number_of_phases_; using Base::saturation_table_number_; - using Base::comp_frac_; using Base::well_index_; using Base::index_of_well_; - using Base::well_type_; using Base::num_components_; using Base::connectionRates_; @@ -398,7 +400,6 @@ namespace Opm void assembleGroupProductionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger); void assembleGroupInjectionControl(const Group& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, const Well::InjectorType& injectorType, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger); - // handle the non reasonable fractions due to numerical overshoot void processFractions() const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index d21c47c9e..d7fb3d92f 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -28,12 +28,16 @@ namespace Opm template StandardWell:: - StandardWell(const Well& well, const int time_step, const Wells* wells, + StandardWell(const Well& well, const int time_step, const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, - const int num_components) - : Base(well, time_step, wells, param, rate_converter, pvtRegionIdx, num_components) + const int num_components, + const int num_phases, + const int index_of_well, + const int first_perf_index, + const std::vector& perf_data) + : Base(well, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, first_perf_index, perf_data) , perf_densities_(number_of_perforations_) , perf_pressure_diffs_(number_of_perforations_) , F0_(numWellConservationEq) @@ -68,7 +72,7 @@ namespace Opm } // counting/updating primary variable numbers - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { // adding a primary variable for water perforation rate per connection numWellEq_ += number_of_perforations_; // adding a primary variable for skin pressure per connection @@ -186,25 +190,33 @@ namespace Opm // Note: currently, the WQTotal definition is still depends on Injector/Producer. assert(comp_idx < num_components_); - if (well_type_ == INJECTOR) { // only single phase injection - // TODO: using comp_frac here is dangerous, it should be changed later - // Most likely, it should be changed to use distr, or at least, we need to update comp_frac_ based on distr - // while solvent might complicate the situation - const auto pu = phaseUsage(); - const int legacyCompIdx = ebosCompIdxToFlowCompIdx(comp_idx); - double comp_frac = 0.0; - if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent - comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent(); - } else if (legacyCompIdx == pu.phase_pos[ Gas ]) { - comp_frac = comp_frac_[legacyCompIdx]; - if (has_solvent) { - comp_frac *= (1.0 - wsolvent()); + if (this->isInjector()) { // only single phase injection + double inj_frac = 0.0; + switch (this->wellEcl().injectorType()) { + case Well::InjectorType::WATER: + if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx))) { + inj_frac = 1.0; } - } else { - comp_frac = comp_frac_[legacyCompIdx]; + break; + case Well::InjectorType::GAS: + if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent + inj_frac = wsolvent(); + } else if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx))) { + inj_frac = has_solvent ? 1.0 - wsolvent() : 1.0; + } + break; + case Well::InjectorType::OIL: + if (comp_idx == int(Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx))) { + inj_frac = 1.0; + } + break; + case Well::InjectorType::MULTI: + // Not supported. + // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + // "Multi phase injectors are not supported, requested for well " + name()); + break; } - - return comp_frac * primary_variables_evaluation_[WQTotal]; + return inj_frac * primary_variables_evaluation_[WQTotal]; } else { // producers return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx); } @@ -364,7 +376,7 @@ namespace Opm const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf]; EvalWell drawdown = pressure - well_pressure; - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; const EvalWell& skin_pressure = primary_variables_evaluation_[pskin_index]; drawdown += skin_pressure; @@ -373,7 +385,7 @@ namespace Opm // producing perforations if ( drawdown.value() > 0 ) { //Do nothing if crossflow is not allowed - if (!allow_cf && well_type_ == INJECTOR) { + if (!allow_cf && this->isInjector()) { return; } @@ -395,7 +407,7 @@ namespace Opm cq_s[oilCompIdx] += vap_oil; // recording the perforation solution gas rate and solution oil rates - if (well_type_ == PRODUCER) { + if (this->isProducer()) { perf_dis_gas_rate = dis_gas.value(); perf_vap_oil_rate = vap_oil.value(); } @@ -403,7 +415,7 @@ namespace Opm } else { //Do nothing if crossflow is not allowed - if (!allow_cf && well_type_ == PRODUCER) { + if (!allow_cf && this->isProducer()) { return; } @@ -471,7 +483,7 @@ namespace Opm } // calculating the perforation solution gas rate and solution oil rates - if (well_type_ == PRODUCER) { + if (this->isProducer()) { if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); @@ -553,12 +565,12 @@ namespace Opm cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); // better way to do here is that use the cq_s and then replace the cq_s_water here? - if (has_polymer && this->has_polymermw && well_type_ == INJECTOR) { + if (has_polymer && this->has_polymermw && this->isInjector()) { handleInjectivityRateAndEquations(intQuants, well_state, perf, cq_s, deferred_logger); } // updating the solution gas rate and solution oil rate - if (well_type_ == PRODUCER) { + if (this->isProducer()) { well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate; well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; } @@ -629,7 +641,7 @@ namespace Opm } // change temperature for injecting fluids - if (well_type_ == INJECTOR && cq_s[activeCompIdx] > 0.0){ + if (this->isInjector() && cq_s[activeCompIdx] > 0.0){ // only handles single phase injection now assert(this->well_ecl_.injectorType() != Well::InjectorType::MULTI); fs.setTemperature(this->well_ecl_.temperature()); @@ -655,7 +667,7 @@ namespace Opm // TODO: the application of well efficiency factor has not been tested with an example yet const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); EvalWell cq_s_poly = cq_s[waterCompIdx] * well_efficiency_factor_; - if (well_type_ == INJECTOR) { + if (this->isInjector()) { cq_s_poly *= wpolymer(); } else { cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection()); @@ -671,7 +683,7 @@ namespace Opm // TODO: the application of well efficiency factor has not been tested with an example yet const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); EvalWell cq_s_foam = cq_s[gasCompIdx] * well_efficiency_factor_; - if (well_type_ == INJECTOR) { + if (this->isInjector()) { cq_s_foam *= wfoam(); } else { cq_s_foam *= extendEval(intQuants.foamConcentration()); @@ -751,7 +763,7 @@ namespace Opm if (wellIsStopped_) { control_eq = getWQTotal(); - } else if (well.isInjector()) { + } else if (this->isInjector()) { const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index]; const auto& controls = well.injectionControls(summaryState); switch(current) { @@ -1309,7 +1321,7 @@ namespace Opm // for injectors, very typical one of the fractions will be one, and it is easy to get zero value // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here - const double relaxation_factor_fractions = (well_type_ == PRODUCER) ? + const double relaxation_factor_fractions = (this->isProducer()) ? relaxationFactorFractionsProducer(old_primary_variables, dwells) : 1.0; @@ -1368,7 +1380,7 @@ namespace Opm updateExtraPrimaryVariables(const BVectorWell& dwells) const { // for the water velocity and skin pressure - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { for (int perf = 0; perf < number_of_perforations_; ++perf) { const int wat_vel_index = Bhp + 1 + perf; const int pskin_index = Bhp + 1 + number_of_perforations_ + perf; @@ -1520,24 +1532,37 @@ namespace Opm // calculate the phase rates based on the primary variables // for producers, this is not a problem, while not sure for injectors here - if (well_type_ == PRODUCER) { + if (this->isProducer()) { const double g_total = primary_variables_[WQTotal]; for (int p = 0; p < number_of_phases_; ++p) { well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = g_total * F[p]; } } else { // injectors - // TODO: using comp_frac_ here is very dangerous, since we do not update it based on the injection phase - // Either we use distr (might conflict with RESV related) or we update comp_frac_ based on the injection phase for (int p = 0; p < number_of_phases_; ++p) { - const double comp_frac = comp_frac_[p]; - well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = comp_frac * primary_variables_[WQTotal]; + well_state.wellRates()[index_of_well_ * number_of_phases_ + p] = 0.0; + } + switch (this->wellEcl().injectorType()) { + case Well::InjectorType::WATER: + well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[Water]] = primary_variables_[WQTotal]; + break; + case Well::InjectorType::GAS: + well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[Gas]] = primary_variables_[WQTotal]; + break; + case Well::InjectorType::OIL: + well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[Oil]] = primary_variables_[WQTotal]; + break; + case Well::InjectorType::MULTI: + // Not supported. + deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + "Multi phase injectors are not supported, requested for well " + name()); + break; } } updateThp(well_state, deferred_logger); // other primary variables related to polymer injectivity study - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { for (int perf = 0; perf < number_of_perforations_; ++perf) { well_state.perfWaterVelocity()[first_perf_ + perf] = primary_variables_[Bhp + 1 + perf]; well_state.perfSkinPressure()[first_perf_ + perf] = primary_variables_[Bhp + 1 + number_of_perforations_ + perf]; @@ -1607,7 +1632,7 @@ namespace Opm return; } - if (well.isInjector() ) + if (this->isInjector() ) { const auto& controls = well.injectionControls(summaryState); @@ -1824,7 +1849,7 @@ namespace Opm // TODO: it only handles the producers for now // the formular for the injectors are not formulated yet - if (well_type_ == INJECTOR) { + if (this->isInjector()) { return; } @@ -1926,7 +1951,7 @@ namespace Opm } // focusing on PRODUCER for now - if (well_type_ == INJECTOR) { + if (this->isInjector()) { return; } @@ -2105,8 +2130,8 @@ namespace Opm // for now, if there is one perforation can produce/inject in the correct // direction, we consider this well can still produce/inject. // TODO: it can be more complicated than this to cause wrong-signed rates - if ( (drawdown < 0. && well_type_ == INJECTOR) || - (drawdown > 0. && well_type_ == PRODUCER) ) { + if ( (drawdown < 0. && this->isInjector()) || + (drawdown > 0. && this->isProducer()) ) { all_drawdown_wrong_direction = false; break; } @@ -2129,15 +2154,15 @@ namespace Opm std::vector well_rates; computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger); - const double sign = (well_type_ == PRODUCER) ? -1. : 1.; + const double sign = (this->isProducer()) ? -1. : 1.; const double threshold = sign * std::numeric_limits::min(); bool can_produce_inject = false; for (const auto value : well_rates) { - if (well_type_ == PRODUCER && value < threshold) { + if (this->isProducer() && value < threshold) { can_produce_inject = true; break; - } else if (well_type_ == INJECTOR && value > threshold) { + } else if (this->isInjector() && value > threshold) { can_produce_inject = true; break; } @@ -2293,7 +2318,6 @@ namespace Opm const std::vector& surf_dens_perf) { // Verify that we have consistent input. - const int np = number_of_phases_; const int nperf = number_of_perforations_; const int num_comp = num_components_; @@ -2336,13 +2360,44 @@ namespace Opm mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate); } } else { + std::fill(mix.begin(), mix.end(), 0.0); // No flow => use well specified fractions for mix. - for (int component = 0; component < num_comp; ++component) { - if (component < np) { - mix[component] = comp_frac_[ ebosCompIdxToFlowCompIdx(component)]; + if (this->isInjector()) { + switch (this->wellEcl().injectorType()) { + case Well::InjectorType::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + case Well::InjectorType::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case Well::InjectorType::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case Well::InjectorType::MULTI: + // Not supported. + // deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + // "Multi phase injectors are not supported, requested for well " + name()); + break; } + } else { + assert(this->isProducer()); + // Using the preferred phase to decide the mix initialization. + switch (this->wellEcl().getPreferredPhase()) { + case Phase::OIL: + mix[FluidSystem::oilCompIdx] = 1.0; + break; + case Phase::GAS: + mix[FluidSystem::gasCompIdx] = 1.0; + break; + case Phase::WATER: + mix[FluidSystem::waterCompIdx] = 1.0; + break; + default: + // No others supported. + break; + } + } - // intialize 0.0 for comIdx >= np; } // Compute volume ratio. x = mix; @@ -2848,6 +2903,7 @@ namespace Opm const int well_index = index_of_well_; const int np = number_of_phases_; + const auto& pu = phaseUsage(); // the weighted total well rate double total_well_rate = 0.0; @@ -2857,12 +2913,22 @@ namespace Opm // Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate // under surface condition is used here - if (well_type_ == INJECTOR) { - primary_variables_[WQTotal] = 0.; - for (int p = 0; p < np; ++p) { - // TODO: the use of comp_frac_ here is dangerous, since the injection phase can be different from - // prefered phasse in WELSPECS, while comp_frac_ only reflect the one specified in WELSPECS - primary_variables_[WQTotal] += well_state.wellRates()[np * well_index + p] * comp_frac_[p]; + if (this->isInjector()) { + switch (this->wellEcl().injectorType()) { + case Well::InjectorType::WATER: + primary_variables_[WQTotal] = well_state.wellRates()[np * well_index + pu.phase_pos[Water]]; + break; + case Well::InjectorType::GAS: + primary_variables_[WQTotal] = well_state.wellRates()[np * well_index + pu.phase_pos[Gas]]; + break; + case Well::InjectorType::OIL: + primary_variables_[WQTotal] = well_state.wellRates()[np * well_index + pu.phase_pos[Oil]]; + break; + case Well::InjectorType::MULTI: + // Not supported. + deferred_logger.warning("MULTI_PHASE_INJECTOR_NOT_SUPPORTED", + "Multi phase injectors are not supported, requested for well " + name()); + break; } } else { for (int p = 0; p < np; ++p) { @@ -2871,7 +2937,6 @@ namespace Opm } if (std::abs(total_well_rate) > 0.) { - const auto pu = phaseUsage(); if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / total_well_rate; } @@ -2882,7 +2947,7 @@ namespace Opm primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / total_well_rate ; } } else { // total_well_rate == 0 - if (well_type_ == INJECTOR) { + if (this->isInjector()) { auto phase = well_ecl_.getInjectionProperties().injectorType; // only single phase injection handled if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { @@ -2907,7 +2972,7 @@ namespace Opm // TODO: it is possible to leave injector as a oil well, // when F_w and F_g both equals to zero, not sure under what kind of circumstance // this will happen. - } else if (well_type_ == PRODUCER) { // producers + } else if (this->isProducer()) { // producers // TODO: the following are not addressed for the solvent case yet if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { primary_variables_[WFrac] = 1.0 / np; @@ -2925,7 +2990,7 @@ namespace Opm primary_variables_[Bhp] = well_state.bhp()[index_of_well_]; // other primary variables related to polymer injection - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { for (int perf = 0; perf < number_of_perforations_; ++perf) { primary_variables_[Bhp + 1 + perf] = well_state.perfWaterVelocity()[first_perf_ + perf]; primary_variables_[Bhp + 1 + number_of_perforations_ + perf] = well_state.perfSkinPressure()[first_perf_ + perf]; @@ -2965,14 +3030,14 @@ namespace Opm // TODO: it is possible it should be a Evaluation const double rho = perf_densities_[0]; - if (well.isInjector() ) + if (this->isInjector() ) { const auto& controls = well.injectionControls(summaryState); const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit) - dp; } - else if (well.isProducer()) { + else if (this->isProducer()) { const auto& controls = well.productionControls(summaryState); const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); @@ -3005,14 +3070,14 @@ namespace Opm const double rho = perf_densities_[0]; double thp = 0.0; - if (well_type_ == INJECTOR) { + if (this->isInjector()) { const int table_id = well_ecl_.vfp_table_number(); const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); } - else if (well_type_ == PRODUCER) { + else if (this->isProducer()) { const int table_id = well_ecl_.vfp_table_number(); const double alq = well_ecl_.alq_value(); const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); @@ -3052,7 +3117,7 @@ namespace Opm // TODO: not sure should based on the well type or injecting/producing peforations // it can be different for crossflow - if (well_type_ == INJECTOR) { + if (this->isInjector()) { // assume fully mixing within injecting wellbore const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex()); const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); @@ -3062,7 +3127,7 @@ namespace Opm if (PolymerModule::hasPlyshlog()) { // we do not calculate the shear effects for injection wells when they do not // inject polymer. - if (well_type_ == INJECTOR && wpolymer() == 0.) { + if (this->isInjector() && wpolymer() == 0.) { return; } // compute the well water velocity with out shear effects. @@ -3398,7 +3463,7 @@ namespace Opm StandardWell:: updateWaterThroughput(const double dt, WellState &well_state) const { - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { for (int perf = 0; perf < number_of_perforations_; ++perf) { const double perf_water_vel = primary_variables_[Bhp + 1 + perf]; // we do not consider the formation damage due to water flowing from reservoir into wellbore @@ -3477,13 +3542,12 @@ namespace Opm using CR = ConvergenceReport; CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid; - const auto& well = well_ecl_; const int well_index = index_of_well_; if (wellIsStopped_) { ctrltype = CR::WellFailure::Type::ControlRate; control_tolerance = 1.e-6; // use smaller tolerance for zero control? } - else if (well.isInjector() ) + else if (this->isInjector() ) { const Opm::Well::InjectorCMode& current = well_state.currentInjectionControls()[well_index]; switch(current) { @@ -3508,7 +3572,7 @@ namespace Opm OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger); } } - else if (well.isProducer() ) + else if (this->isProducer() ) { const Well::ProducerCMode& current = well_state.currentProductionControls()[well_index]; switch(current) { @@ -3563,7 +3627,7 @@ namespace Opm // if different types of extra equations are involved, this function needs to be refactored further // checking the convergence of the extra equations related to polymer injectivity - if (this->has_polymermw && well_type_ == INJECTOR) { + if (this->has_polymermw && this->isInjector()) { // checking the convergence of the perforation rates const double wat_vel_tol = 1.e-8; const int dummy_component = -1; @@ -3612,7 +3676,7 @@ namespace Opm { // the source term related to transport of molecular weight EvalWell cq_s_polymw = cq_s_poly; - if (well_type_ == INJECTOR) { + if (this->isInjector()) { const int wat_vel_index = Bhp + 1 + perf; const EvalWell water_velocity = primary_variables_evaluation_[wat_vel_index]; if (water_velocity > 0.) { // injecting @@ -3624,7 +3688,7 @@ namespace Opm // going-back to the wellbore through injector cq_s_polymw *= 0.; } - } else if (well_type_ == PRODUCER) { + } else if (this->isProducer()) { if (cq_s_polymw < 0.) { cq_s_polymw *= extendEval(int_quants.polymerMoleWeight() ); } else { diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index ea5cb65e2..541d68c30 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -118,11 +118,15 @@ namespace Opm compositionSwitchEnabled, Indices::numPhases >; /// Constructor - WellInterface(const Well& well, const int time_step, const Wells* wells, + WellInterface(const Well& well, const int time_step, const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, - const int num_components); + const int num_components, + const int num_phases, + const int index_of_well, + const int first_perf_index, + const std::vector& perf_data); /// Virutal destructor virtual ~WellInterface() {} @@ -130,15 +134,18 @@ namespace Opm /// Well name. const std::string& name() const; + /// True if the well is an injector. + bool isInjector() const; + + /// True if the well is a producer. + bool isProducer() const; + /// Index of well in the wells struct and wellState int indexOfWell() const; /// Well cells. const std::vector& cells() const {return well_cells_; } - /// Well type, INJECTOR or PRODUCER. - WellType wellType() const; - /// Well controls WellControls* wellControls() const; @@ -286,30 +293,12 @@ namespace Opm const int current_step_; - // the index of well in Wells struct - int index_of_well_; - // simulation parameters const ModelParameters& param_; - // well type - // INJECTOR or PRODUCER - enum WellType well_type_; - - // number of phases - int number_of_phases_; - - // component fractions for each well - // typically, it should apply to injection wells - std::vector comp_frac_; - // number of the perforations for this well int number_of_perforations_; - // record the index of the first perforation - // of states of individual well. - int first_perf_; - // well index for each perforation std::vector well_index_; @@ -372,6 +361,18 @@ namespace Opm const int num_components_; + // number of phases + int number_of_phases_; + + // the index of well in Wells struct + int index_of_well_; + + // record the index of the first perforation + // of states of individual well. + int first_perf_; + + std::vector originalConnectionIndex_; + std::vector connectionRates_; bool wellIsStopped_; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index a3bd8d6d2..cdb449b52 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -27,75 +27,58 @@ namespace Opm template WellInterface:: - WellInterface(const Well& well, const int time_step, const Wells* wells, + WellInterface(const Well& well, + const int time_step, const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, - const int num_components) + const int num_components, + const int num_phases, + const int index_of_well, + const int first_perf_index, + const std::vector& perf_data) : well_ecl_(well) , current_step_(time_step) , param_(param) , rateConverter_(rate_converter) , pvtRegionIdx_(pvtRegionIdx) , num_components_(num_components) + , number_of_phases_(num_phases) + , index_of_well_(index_of_well) + , first_perf_(first_perf_index) { if (time_step < 0) { OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface"); } - if (!wells) { - OPM_THROW(std::invalid_argument, "Null pointer of Wells is used to construct WellInterface"); - } + ref_depth_ = well.getRefDepth(); - const std::string& well_name = well.name(); - - // looking for the location of the well in the wells struct - int index_well; - for (index_well = 0; index_well < wells->number_of_wells; ++index_well) { - if (well_name == std::string(wells->name[index_well])) { - break; - } - } - - // should not enter the constructor if the well does not exist in the wells struct - // here, just another assertion. - assert(index_well != wells->number_of_wells); - - index_of_well_ = index_well; - well_type_ = wells->type[index_well]; - number_of_phases_ = wells->number_of_phases; - - // copying the comp_frac - { - comp_frac_.resize(number_of_phases_); - const int index_begin = index_well * number_of_phases_; - std::copy(wells->comp_frac + index_begin, - wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() ); - } - - ref_depth_ = wells->depth_ref[index_well]; + // We do not want to count SHUT perforations here, so + // it would be wrong to use wells.getConnections().size(). + number_of_perforations_ = perf_data.size(); // perforations related { - const int perf_index_begin = wells->well_connpos[index_well]; - const int perf_index_end = wells->well_connpos[index_well + 1]; - number_of_perforations_ = perf_index_end - perf_index_begin; - first_perf_ = perf_index_begin; - well_cells_.resize(number_of_perforations_); - std::copy(wells->well_cells + perf_index_begin, - wells->well_cells + perf_index_end, - well_cells_.begin() ); - well_index_.resize(number_of_perforations_); - std::copy(wells->WI + perf_index_begin, - wells->WI + perf_index_end, - well_index_.begin() ); - saturation_table_number_.resize(number_of_perforations_); - std::copy(wells->sat_table_id + perf_index_begin, - wells->sat_table_id + perf_index_end, - saturation_table_number_.begin() ); + int perf = 0; + for (const auto& pd : perf_data) { + well_cells_[perf] = pd.cell_index; + well_index_[perf] = pd.connection_transmissibility_factor; + saturation_table_number_[perf] = pd.satnum_id; + ++perf; + } + + int all_perf = 0; + originalConnectionIndex_.reserve(perf_data.size()); + for (const auto connection : well.getConnections()) { + if (connection.state() == Connection::State::OPEN) { + originalConnectionIndex_.push_back(all_perf); + } + ++all_perf; + } + assert(originalConnectionIndex_.size() == perf_data.size()); } // initialization of the completions mapping @@ -205,11 +188,22 @@ namespace Opm template - WellType + bool WellInterface:: - wellType() const + isInjector() const { - return well_type_; + return well_ecl_.isInjector(); + } + + + + + template + bool + WellInterface:: + isProducer() const + { + return well_ecl_.isProducer(); } @@ -781,7 +775,7 @@ namespace Opm { // currently, we only updateWellTestState for producers - if (wellType() != PRODUCER) { + if (this->isInjector()) { return; } @@ -930,7 +924,8 @@ namespace Opm bool allCompletionsClosed = true; const auto& connections = well_ecl_.getConnections(); for (const auto& connection : connections) { - if (!well_test_state.hasCompletion(name(), connection.complnum())) { + if (connection.state() == Connection::State::OPEN + && !well_test_state.hasCompletion(name(), connection.complnum())) { allCompletionsClosed = false; } } @@ -1166,7 +1161,7 @@ namespace Opm // we need to get the table number through the parser, in case THP constraint/target is not there. // When THP control/limit is not active, if available VFP table is provided, we will still need to // update THP value. However, it will only used for output purpose. - if (well_type_ == PRODUCER) { // producer + if (isProducer()) { // producer const int table_id = well_ecl_.vfp_table_number(); if (table_id <= 0) { return false; @@ -1264,10 +1259,12 @@ namespace Opm const auto& connections = well_ecl_.getConnections(); int perfIdx = 0; for (const auto& connection : connections) { - if (wellTestState.hasCompletion(name(), connection.complnum())) { - well_index_[perfIdx] = 0.0; + if (connection.state() == Connection::State::OPEN) { + if (wellTestState.hasCompletion(name(), connection.complnum())) { + well_index_[perfIdx] = 0.0; + } + perfIdx++; } - perfIdx++; } } @@ -1295,7 +1292,7 @@ namespace Opm void WellInterface::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) { - const auto& connection = well_ecl_.getConnections()[perfIdx]; + const auto& connection = well_ecl_.getConnections()[originalConnectionIndex_[perfIdx]]; if (well_ecl_.getDrainageRadius() < 0) { if (new_well && perfIdx == 0) { deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero"); diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 03673615c..f963d910c 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -24,6 +24,8 @@ #include #include #include +#include +#include #include #include @@ -42,139 +44,54 @@ namespace Opm typedef std::array< int, 3 > mapentry_t; typedef std::map< std::string, mapentry_t > WellMapType; - template - void init(const Wells* wells, const State& state) - { - init(wells, state.pressure()); - } - /// Allocate and initialize if wells is non-null. /// Also tries to give useful initial values to the bhp() and /// wellRates() fields, depending on controls. The /// perfRates() field is filled with zero, and perfPress() /// with -1e100. - void init(const Wells* wells, const std::vector& cellPressures) + void init(const std::vector& cellPressures, + const std::vector& wells_ecl, + const PhaseUsage& pu, + const std::vector>& well_perf_data, + const SummaryState& summary_state) { // clear old name mapping wellMap_.clear(); - wells_.reset( clone_wells( wells ) ); - if (wells) { - const int nw = wells->number_of_wells; - const int np = wells->number_of_phases; - bhp_.resize(nw); - thp_.resize(nw); + well_perf_data_ = well_perf_data; + + { + // const int nw = wells->number_of_wells; + const int nw = wells_ecl.size(); + // const int np = wells->number_of_phases; + const int np = pu.num_phases; + bhp_.resize(nw, 0.0); + thp_.resize(nw, 0.0); temperature_.resize(nw, 273.15 + 20); // standard temperature for now wellrates_.resize(nw * np, 0.0); + int connpos = 0; for (int w = 0; w < nw; ++w) { - assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER)); - const WellControls* ctrl = wells->ctrls[w]; - const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w]; + const Well& well = wells_ecl[w]; - // setup wellname -> well index mapping - { - assert( wells->name[ w ] ); - std::string name( wells->name[ w ] ); - assert( name.size() > 0 ); - mapentry_t& wellMapEntry = wellMap_[name]; - wellMapEntry[ 0 ] = w; - wellMapEntry[ 1 ] = wells->well_connpos[w]; - // also store the number of perforations in this well - wellMapEntry[ 2 ] = num_perf_this_well; - } + // Initialize bhp(), thp(), wellRates(). + initSingleWell(cellPressures, w, well, pu, summary_state); - if ( num_perf_this_well == 0 ) - { - // No perforations of the well. Initialize to zero. - for (int p = 0; p < np; ++p) { - wellrates_[np*w + p] = 0.0; - } - bhp_[w] = 0.; - thp_[w] = 0.; - continue; - } - - if (well_controls_well_is_stopped(ctrl)) { - // Stopped well: - // 1. Rates: assign zero well rates. - for (int p = 0; p < np; ++p) { - wellrates_[np*w + p] = 0.0; - } - // 2. Bhp: assign bhp equal to bhp control, if - // applicable, otherwise assign equal to - // first perforation cell pressure. - if (well_controls_get_current_type(ctrl) == BHP) { - bhp_[w] = well_controls_get_current_target( ctrl ); - } else { - const int first_cell = wells->well_cells[wells->well_connpos[w]]; - bhp_[w] = cellPressures[first_cell]; - } - } else if (well_controls_get_current(ctrl) == -1) { - // Well under group control. - // 1. Rates: assign zero well rates. - for (int p = 0; p < np; ++p) { - wellrates_[np*w + p] = 0.0; - } - // 2. Bhp: initialize bhp to be a - // little above or below (depending on if - // the well is an injector or producer) - // pressure in first perforation cell. - const int first_cell = wells->well_cells[wells->well_connpos[w]]; - const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99; - bhp_[w] = safety_factor*cellPressures[first_cell]; - } else { - // Open well: - // 1. Rates: initialize well rates to match controls - // if type is SURFACE_RATE. Otherwise, we - // cannot set the correct value here, so we - // assign a small rate with the correct - // sign so that any logic depending on that - // sign will work as expected. - if (well_controls_get_current_type(ctrl) == SURFACE_RATE) { - const double rate_target = well_controls_get_current_target(ctrl); - const double * distr = well_controls_get_current_distr( ctrl ); - for (int p = 0; p < np; ++p) { - wellrates_[np*w + p] = rate_target * distr[p]; - } - } else { - const double small_rate = 0.0; //1e-14; - const double sign = (wells->type[w] == INJECTOR) ? 1.0 : -1.0; - for (int p = 0; p < np; ++p) { - wellrates_[np*w + p] = small_rate * sign; - } - } - - // 2. Bhp: initialize bhp to be target pressure if - // bhp-controlled well, otherwise set to a - // little above or below (depending on if - // the well is an injector or producer) - // pressure in first perforation cell. - if (well_controls_get_current_type(ctrl) == BHP) { - bhp_[w] = well_controls_get_current_target( ctrl ); - } else { - const int first_cell = wells->well_cells[wells->well_connpos[w]]; - const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99; - bhp_[w] = safety_factor*cellPressures[first_cell]; - } - } - - // 3. Thp: assign thp equal to thp target/limit, if applicable, - // otherwise keep it zero. Basically, the value should not be used - // in the simulation at all. - const int nwc = well_controls_get_num(ctrl); - for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(ctrl, ctrl_index) == THP) { - thp_[w] = well_controls_iget_target(ctrl, ctrl_index); - break; - } - } + // Setup wellname -> well index mapping. + const int num_perf_this_well = well_perf_data[w].size(); + std::string name = well.name(); + assert( name.size() > 0 ); + mapentry_t& wellMapEntry = wellMap_[name]; + wellMapEntry[ 0 ] = w; + wellMapEntry[ 1 ] = connpos; + wellMapEntry[ 2 ] = num_perf_this_well; + connpos += num_perf_this_well; } // The perforation rates and perforation pressures are // not expected to be consistent with bhp_ and wellrates_ // after init(). - perfrates_.resize(wells->well_connpos[nw], 0.0); - perfpress_.resize(wells->well_connpos[nw], -1e100); + perfrates_.resize(connpos, 0.0); + perfpress_.resize(connpos, -1e100); } } @@ -264,14 +181,11 @@ namespace Opm well.rates.set( rt::gas, wv[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ] ); } - const int num_perf_well = this->wells_->well_connpos[ well_index + 1 ] - - this->wells_->well_connpos[ well_index ]; + const int num_perf_well = this->well_perf_data_[well_index].size(); well.connections.resize(num_perf_well); for( int i = 0; i < num_perf_well; ++i ) { - const auto wi = this->wells_->well_connpos[ well_index ] + i; - const auto active_index = this->wells_->well_cells[ wi ]; - + const auto active_index = this->well_perf_data_[well_index][i].cell_index; auto& connection = well.connections[ i ]; connection.index = globalCellIdxMap[active_index]; connection.pressure = this->perfPress()[ itr.second[1] + i ]; @@ -284,32 +198,10 @@ namespace Opm } - virtual ~WellState() {} - + virtual ~WellState() = default; WellState() = default; - WellState( const WellState& rhs ) : - bhp_( rhs.bhp_ ), - thp_( rhs.thp_ ), - temperature_( rhs.temperature_ ), - wellrates_( rhs.wellrates_ ), - perfrates_( rhs.perfrates_ ), - perfpress_( rhs.perfpress_ ), - wellMap_( rhs.wellMap_ ), - wells_( clone_wells( rhs.wells_.get() ) ) - {} - - WellState& operator=( const WellState& rhs ) { - this->bhp_ = rhs.bhp_; - this->thp_ = rhs.thp_; - this->temperature_ = rhs.temperature_; - this->wellrates_ = rhs.wellrates_; - this->perfrates_ = rhs.perfrates_; - this->perfpress_ = rhs.perfpress_; - this->wellMap_ = rhs.wellMap_; - this->wells_.reset( clone_wells( rhs.wells_.get() ) ); - - return *this; - } + WellState(const WellState& rhs) = default; + WellState& operator=(const WellState& rhs) = default; private: std::vector bhp_; @@ -321,11 +213,139 @@ namespace Opm WellMapType wellMap_; + void initSingleWell(const std::vector& cellPressures, + const int w, + const Well& well, + const PhaseUsage& pu, + const SummaryState& summary_state) + { + assert(well.isInjector() || well.isProducer()); + + // Set default zero initial well rates. + // May be overwritten below. + const int np = pu.num_phases; + for (int p = 0; p < np; ++p) { + wellrates_[np*w + p] = 0.0; + } + + const int num_perf_this_well = well_perf_data_[w].size(); + if ( num_perf_this_well == 0 ) { + // No perforations of the well. Initialize to zero. + bhp_[w] = 0.; + thp_[w] = 0.; + return; + } + + const auto inj_controls = well.isInjector() ? well.injectionControls(summary_state) : Well::InjectionControls(0); + const auto prod_controls = well.isProducer() ? well.productionControls(summary_state) : Well::ProductionControls(0); + + const bool is_bhp = well.isInjector() ? (inj_controls.cmode == Well::InjectorCMode::BHP) + : (prod_controls.cmode == Well::ProducerCMode::BHP); + const double bhp_limit = well.isInjector() ? inj_controls.bhp_limit : prod_controls.bhp_limit; + const bool is_grup = well.isInjector() ? (inj_controls.cmode == Well::InjectorCMode::GRUP) + : (prod_controls.cmode == Well::ProducerCMode::GRUP); + + const double inj_surf_rate = well.isInjector() ? inj_controls.surface_rate : 0.0; // To avoid a "maybe-uninitialized" warning. + + if (well.getStatus() == Well::Status::STOP) { + // Stopped well: + // 1. Rates: zero well rates. + // 2. Bhp: assign bhp equal to bhp control, if + // applicable, otherwise assign equal to + // first perforation cell pressure. + if (is_bhp) { + bhp_[w] = bhp_limit; + } else { + const int first_cell = well_perf_data_[w][0].cell_index; + bhp_[w] = cellPressures[first_cell]; + } + } else if (is_grup) { + // Well under group control. + // 1. Rates: zero well rates. + // 2. Bhp: initialize bhp to be a + // little above or below (depending on if + // the well is an injector or producer) + // pressure in first perforation cell. + const int first_cell = well_perf_data_[w][0].cell_index; + const double safety_factor = well.isInjector() ? 1.01 : 0.99; + bhp_[w] = safety_factor*cellPressures[first_cell]; + } else { + // Open well, under own control: + // 1. Rates: initialize well rates to match + // controls if type is ORAT/GRAT/WRAT + // (producer) or RATE (injector). + // Otherwise, we cannot set the correct + // value here and initialize to zero rate. + if (well.isInjector()) { + if (inj_controls.cmode == Well::InjectorCMode::RATE) { + switch (inj_controls.injector_type) { + case Well::InjectorType::WATER: + assert(pu.phase_used[BlackoilPhases::Aqua]); + wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate; + break; + case Well::InjectorType::GAS: + assert(pu.phase_used[BlackoilPhases::Vapour]); + wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate; + break; + case Well::InjectorType::OIL: + assert(pu.phase_used[BlackoilPhases::Liquid]); + wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate; + break; + case Well::InjectorType::MULTI: + // Not currently handled, keep zero init. + break; + } + } else { + // Keep zero init. + } + } else { + assert(well.isProducer()); + // Note negative rates for producing wells. + switch (prod_controls.cmode) { + case Well::ProducerCMode::ORAT: + assert(pu.phase_used[BlackoilPhases::Liquid]); + wellrates_[np*w + pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate; + break; + case Well::ProducerCMode::WRAT: + assert(pu.phase_used[BlackoilPhases::Aqua]); + wellrates_[np*w + pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate; + break; + case Well::ProducerCMode::GRAT: + assert(pu.phase_used[BlackoilPhases::Vapour]); + wellrates_[np*w + pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate; + break; + default: + // Keep zero init. + break; + } + } + // 2. Bhp: initialize bhp to be target pressure if + // bhp-controlled well, otherwise set to a + // little above or below (depending on if + // the well is an injector or producer) + // pressure in first perforation cell. + if (is_bhp) { + bhp_[w] = bhp_limit; + } else { + const int first_cell = well_perf_data_[w][0].cell_index; + const double safety_factor = well.isInjector() ? 1.01 : 0.99; + bhp_[w] = safety_factor*cellPressures[first_cell]; + } + } + + // 3. Thp: assign thp equal to thp target/limit, if such a limit exists, + // otherwise keep it zero. + const bool has_thp = well.isInjector() ? inj_controls.hasControl(Well::InjectorCMode::THP) + : prod_controls.hasControl(Well::ProducerCMode::THP); + const double thp_limit = well.isInjector() ? inj_controls.thp_limit : prod_controls.thp_limit; + if (has_thp) { + thp_[w] = thp_limit; + } + + } + protected: - struct wdel { - void operator()( Wells* w ) { destroy_wells( w ); } - }; - std::unique_ptr< Wells, wdel > wells_; + std::vector> well_perf_data_; }; } // namespace Opm diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 70bffdb92..643c7b8ab 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -52,6 +52,8 @@ namespace Opm public: typedef BaseType :: WellMapType WellMapType; + virtual ~WellStateFullyImplicitBlackoil() = default; + // TODO: same definition with WellInterface, eventually they should go to a common header file. static const int Water = BlackoilPhases::Aqua; static const int Oil = BlackoilPhases::Liquid; @@ -67,29 +69,29 @@ namespace Opm /// Allocate and initialize if wells is non-null. Also tries /// to give useful initial values to the bhp(), wellRates() /// and perfPhaseRates() fields, depending on controls - void init(const Wells* wells, - const std::vector& cellPressures, + void init(const std::vector& cellPressures, const Schedule& schedule, const std::vector& wells_ecl, const int report_step, const WellStateFullyImplicitBlackoil* prevState, - const PhaseUsage& pu) + const PhaseUsage& pu, + const std::vector>& well_perf_data, + const SummaryState& summary_state) { // call init on base class - BaseType :: init(wells, cellPressures); + BaseType :: init(cellPressures, wells_ecl, pu, well_perf_data, summary_state); - // if there are no well, do nothing in init - if (wells == 0) { - return; - } - - const int nw = wells->number_of_wells; + const int nw = wells_ecl.size(); if( nw == 0 ) return ; // Initialize perfphaserates_, which must be done here. - const int np = wells->number_of_phases; - const int nperf = wells->well_connpos[nw]; + const int np = pu.num_phases; + + int nperf = 0; + for (const auto& wpd : well_perf_data) { + nperf += wpd.size(); + } well_reservoir_rates_.resize(nw * np, 0.0); well_dissolved_gas_rates_.resize(nw, 0.0); @@ -107,23 +109,9 @@ namespace Opm const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE + ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::INJECTION_UPDATE; - - for (int w = 0; w name[w]); - for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) { - if (well_name == wells_ecl[index_well_ecl].name()) { - break; - } - } - - // It should be able to find in wells_ecl. - if (index_well_ecl == nw_wells_ecl) { - OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl "); - } - - effective_events_occurred_[w] = (schedule.hasWellEvent(well_name, effective_events_mask, report_step) ); + for (int w = 0; w < nw; ++w) { + effective_events_occurred_[w] + = schedule.hasWellEvent(wells_ecl[w].name(), effective_events_mask, report_step); } } // end of if (!well_ecl.empty() ) @@ -139,21 +127,20 @@ namespace Opm perf_skin_pressure_.clear(); perf_skin_pressure_.resize(nperf, 0.0); + int connpos = 0; for (int w = 0; w < nw; ++w) { - assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER)); - const WellControls* ctrl = wells->ctrls[w]; - - const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w]; - // Open well: Initialize perfphaserates_ to well + // Initialize perfphaserates_ to well // rates divided by the number of perforations. - for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { - if (well_controls_well_is_open(ctrl)) { + const int num_perf_this_well = well_perf_data[w].size(); + for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) { + if (wells_ecl[w].getStatus() == Well::Status::OPEN) { for (int p = 0; p < np; ++p) { perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well); } } - perfPress()[perf] = cellPressures[wells->well_cells[perf]]; + perfPress()[perf] = cellPressures[well_perf_data[w][perf-connpos].cell_index]; } + connpos += num_perf_this_well; } current_injection_controls_.resize(nw); @@ -166,13 +153,14 @@ namespace Opm // intialize wells that have been there before // order may change so the mapping is based on the well name - if(prevState && !prevState->wellMap().empty()) { - typedef typename WellMapType :: const_iterator const_iterator; - const_iterator end = prevState->wellMap().end(); + if (prevState && !prevState->wellMap().empty()) { + connpos = 0; + auto end = prevState->wellMap().end(); for (int w = 0; w < nw; ++w) { - const std::string name( wells->name[ w ] ); - const_iterator it = prevState->wellMap().find( name ); - if( it != end ) + const Well& well = wells_ecl[w]; + const int num_perf_this_well = well_perf_data[w].size(); + auto it = prevState->wellMap().find(well.name()); + if ( it != end ) { const int oldIndex = (*it).second[ 0 ]; const int newIndex = w; @@ -210,19 +198,18 @@ namespace Opm // perfPhaseRates const int oldPerf_idx_beg = (*it).second[ 1 ]; const int num_perf_old_well = (*it).second[ 2 ]; - const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex]; // copy perforation rates when the number of perforations is equal, // otherwise initialize perfphaserates to well rates divided by the number of perforations. if( num_perf_old_well == num_perf_this_well ) { int old_perf_phase_idx = oldPerf_idx_beg *np; - for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np; - perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx ) + for (int perf_phase_idx = connpos*np; + perf_phase_idx < (connpos + num_perf_this_well)*np; ++perf_phase_idx, ++old_perf_phase_idx ) { perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ]; } } else { - for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) { + for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) { for (int p = 0; p < np; ++p) { perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well); } @@ -232,8 +219,7 @@ namespace Opm if( num_perf_old_well == num_perf_this_well ) { int oldPerf_idx = oldPerf_idx_beg; - for (int perf = wells->well_connpos[ newIndex ]; - perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx ) + for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx ) { perfPress()[ perf ] = prevState->perfPress()[ oldPerf_idx ]; } @@ -243,8 +229,7 @@ namespace Opm if( num_perf_old_well == num_perf_this_well ) { int oldPerf_idx = oldPerf_idx_beg; - for (int perf = wells->well_connpos[ newIndex ]; - perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx ) + for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx ) { perfRateSolvent()[ perf ] = prevState->perfRateSolvent()[ oldPerf_idx ]; } @@ -258,8 +243,7 @@ namespace Opm if( num_perf_old_well == num_perf_this_well ) { int oldPerf_idx = oldPerf_idx_beg; - for (int perf = wells->well_connpos[ newIndex ]; - perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx ) + for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf, ++oldPerf_idx ) { perf_water_throughput_[ perf ] = prevState->perfThroughput()[ oldPerf_idx ]; perf_skin_pressure_[ perf ] = prevState->perfSkinPressure()[ oldPerf_idx ]; @@ -269,21 +253,16 @@ namespace Opm } } - // If in the new step, there is no THP related target/limit anymore, its thp value should be // set to zero. - const WellControls* ctrl = wells->ctrls[w]; - const int nwc = well_controls_get_num(ctrl); - int ctrl_index = 0; - for (; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(ctrl, ctrl_index) == THP) { - break; - } - } - // not finding any thp related control/limits - if (ctrl_index == nwc) { - thp()[w] = 0.; + const bool has_thp = well.isInjector() ? well.injectionControls(summary_state).hasControl(Well::InjectorCMode::THP) + : well.productionControls(summary_state).hasControl(Well::ProducerCMode::THP); + if (!has_thp) { + thp()[w] = 0.0; } + + // Increment connection position offset. + connpos += num_perf_this_well; } } @@ -303,14 +282,19 @@ namespace Opm } - void resize(const Wells* wells, const std::vector& wells_ecl, const Schedule& schedule, - const bool handle_ms_well, const size_t numCells, const PhaseUsage& pu) + void resize(const std::vector& wells_ecl, + const Schedule& schedule, + const bool handle_ms_well, + const size_t numCells, + const PhaseUsage& pu, + const std::vector>& well_perf_data, + const SummaryState& summary_state) { const std::vector tmp(numCells, 0.0); // <- UGLY HACK to pass the size - init(wells, tmp, schedule, wells_ecl, 0, nullptr, pu); + init(tmp, schedule, wells_ecl, 0, nullptr, pu, well_perf_data, summary_state); if (handle_ms_well) { - initWellStateMSWell(wells, wells_ecl, pu, nullptr); + initWellStateMSWell(wells_ecl, pu, nullptr); } } @@ -434,13 +418,6 @@ namespace Opm phs.at( pu.phase_pos[Gas] ) = rt::gas; } - if (pu.has_solvent) { - // add solvent component - for( int w = 0; w < nw; ++w ) { - res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) ); - } - } - /* this is a reference or example on **how** to convert from * WellState to something understood by opm-output. it is intended * to be properly implemented and maintained as a part of @@ -491,10 +468,14 @@ namespace Opm well.rates.set( rt::well_potential_gas, this->well_potentials_[well_rate_index + pu.phase_pos[Gas]] ); } + if ( pu.has_solvent ) { + well.rates.set( rt::solvent, solventWellRate(w) ); + } + well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] ); well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] ); - int local_comp_index = 0; + size_t local_comp_index = 0; for( auto& comp : well.connections) { const auto rates = this->perfPhaseRates().begin() + (np * wt.second[ 1 ]) @@ -505,7 +486,7 @@ namespace Opm comp.rates.set( phs[ i ], *(rates + i) ); } } - assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]); + assert(local_comp_index == this->well_perf_data_[w].size()); const auto nseg = this->numSegments(w); for (auto seg_ix = 0*nseg; seg_ix < nseg; ++seg_ix) { @@ -520,11 +501,11 @@ namespace Opm /// init the MS well related. - void initWellStateMSWell(const Wells* wells, const std::vector& wells_ecl, + void initWellStateMSWell(const std::vector& wells_ecl, const PhaseUsage& pu, const WellStateFullyImplicitBlackoil* prev_well_state) { // still using the order in wells - const int nw = wells->number_of_wells; + const int nw = wells_ecl.size(); if (nw == 0) { return; } @@ -538,24 +519,12 @@ namespace Opm seg_number_.clear(); nseg_ = 0; + int connpos = 0; // in the init function, the well rates and perforation rates have been initialized or copied from prevState // what we do here, is to set the segment rates and perforation rates for (int w = 0; w < nw; ++w) { - const int nw_wells_ecl = wells_ecl.size(); - int index_well_ecl = 0; - const std::string well_name(wells->name[w]); - for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) { - if (well_name == wells_ecl[index_well_ecl].name()) { - break; - } - } - - // It should be able to find in wells_ecl. - if (index_well_ecl == nw_wells_ecl) { - OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl "); - } - - const auto& well_ecl = wells_ecl[index_well_ecl]; + const auto& well_ecl = wells_ecl[w]; + int num_perf_this_well = well_perf_data_[w].size(); top_segment_index_.push_back(nseg_); if ( !well_ecl.isMultiSegment() ) { // not multi-segment well nseg_ += 1; @@ -571,7 +540,6 @@ namespace Opm const WellConnections& completion_set = well_ecl.getConnections(); // number of segment for this single well const int well_nseg = segment_set.size(); - // const int nperf = completion_set.size(); int n_activeperf = 0; nseg_ += well_nseg; for (auto segID = 0*well_nseg; segID < well_nseg; ++segID) { @@ -605,8 +573,8 @@ namespace Opm // for the segrates_, now it becomes a recursive solution procedure. { const int np = numPhases(); - const int start_perf = wells->well_connpos[w]; - const int start_perf_next_well = wells->well_connpos[w + 1]; + const int start_perf = connpos; + const int start_perf_next_well = connpos + num_perf_this_well; // make sure the information from wells_ecl consistent with wells assert(n_activeperf == (start_perf_next_well - start_perf)); @@ -640,7 +608,7 @@ namespace Opm // top segment is always the first one, and its pressure is the well bhp segpress_.push_back(bhp()[w]); const int top_segment = top_segment_index_[w]; - const int start_perf = wells->well_connpos[w]; + const int start_perf = connpos; for (int seg = 1; seg < well_nseg; ++seg) { if ( !segment_perforations[seg].empty() ) { const int first_perf = segment_perforations[seg][0]; @@ -654,6 +622,7 @@ namespace Opm } } } + connpos += num_perf_this_well; } assert(int(segpress_.size()) == nseg_); assert(int(segrates_.size()) == nseg_ * numPhases() ); @@ -663,8 +632,7 @@ namespace Opm const auto& end = prev_well_state->wellMap().end(); const int np = numPhases(); for (int w = 0; w < nw; ++w) { - const std::string name( wells->name[w] ); - const auto& it = prev_well_state->wellMap().find( name ); + const auto& it = prev_well_state->wellMap().find( wells_ecl[w].name() ); if (it != end) { // the well is found in the prev_well_state // TODO: the well with same name can change a lot, like they might not have same number of segments @@ -736,8 +704,13 @@ namespace Opm /// One rate pr well double solventWellRate(const int w) const { + int connpos = 0; + for (int iw = 0; iw < w; ++iw) { + connpos += this->well_perf_data_[iw].size(); + } double solvent_well_rate = 0.0; - for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) { + const int endperf = connpos + this->well_perf_data_[w].size(); + for (int perf = connpos; perf < endperf; ++perf ) { solvent_well_rate += perfRateSolvent_[perf]; } return solvent_well_rate; diff --git a/tests/test_stoppedwells.cpp b/tests/test_stoppedwells.cpp index 804c71e30..804c1de6a 100644 --- a/tests/test_stoppedwells.cpp +++ b/tests/test_stoppedwells.cpp @@ -27,12 +27,6 @@ #include #include -#include -#include -#include -#include -#include - using namespace Opm; @@ -42,54 +36,24 @@ BOOST_AUTO_TEST_CASE(TestStoppedWells) Opm::Parser parser; Opm::Deck deck(parser.parseFile(filename)); Opm::EclipseState eclipseState(deck); - Opm::GridManager vanguard(eclipseState.getInputGrid()); const auto& grid = eclipseState.getInputGrid(); const TableManager table ( deck ); const Eclipse3DProperties eclipseProperties ( deck , table, grid); const Opm::Runspec runspec (deck); const Schedule sched(deck, grid, eclipseProperties, runspec); - Opm::SummaryState summaryState(std::chrono::system_clock::from_time_t(sched.getStartTime())); - - double target_surfacerate_inj; - double target_surfacerate_prod; - - const std::vector pressure = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; // Both wells are open in the first schedule step { - Opm::WellsManager wellsManager(eclipseState , sched, summaryState, 0, *vanguard.c_grid()); - const Wells* wells = wellsManager.c_wells(); - const struct WellControls* ctrls0 = wells->ctrls[0]; - const struct WellControls* ctrls1 = wells->ctrls[1]; - BOOST_CHECK(well_controls_well_is_open(ctrls0)); - BOOST_CHECK(well_controls_well_is_open(ctrls1)); - - target_surfacerate_inj = well_controls_iget_target(ctrls0 , 0); - target_surfacerate_prod = well_controls_iget_target(ctrls1 , 0); - - WellState wellstate; - wellstate.init(wells, pressure); - const std::vector wellrates = wellstate.wellRates(); - BOOST_CHECK_EQUAL (target_surfacerate_inj, wellrates[2]); // Gas injector - BOOST_CHECK_EQUAL (target_surfacerate_prod, wellrates[4]); // Oil target rate + auto wells = sched.getWells(0); + BOOST_CHECK(wells[0].getStatus() == Opm::Well::Status::OPEN); + BOOST_CHECK(wells[1].getStatus() == Opm::Well::Status::OPEN); } // The injector is stopped { - Opm::WellsManager wellsManager(eclipseState, sched, summaryState, 1 , *vanguard.c_grid()); - const Wells* wells = wellsManager.c_wells(); - const struct WellControls* ctrls0 = wells->ctrls[0]; - const struct WellControls* ctrls1 = wells->ctrls[1]; - BOOST_CHECK(well_controls_well_is_stopped(ctrls0)); // injector is stopped - BOOST_CHECK(well_controls_well_is_open(ctrls1)); - - WellState wellstate; - wellstate.init(wells, pressure); - - const std::vector wellrates = wellstate.wellRates(); - BOOST_CHECK_EQUAL (0, wellrates[2]); // Gas injector - BOOST_CHECK_EQUAL (target_surfacerate_prod, wellrates[4]); // Oil target rate + auto wells = sched.getWells(1); + BOOST_CHECK(wells[0].getStatus() == Opm::Well::Status::STOP); + BOOST_CHECK(wells[1].getStatus() == Opm::Well::Status::OPEN); } - } diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index f81761b9d..d1657ca77 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -77,36 +77,13 @@ struct SetupTest { schedule.reset( new Opm::Schedule(deck, ecl_state->getInputGrid(), eclipseProperties, runspec)); summaryState.reset( new Opm::SummaryState(std::chrono::system_clock::from_time_t(schedule->getStartTime()))); } - - // Create grid. - const std::vector& porv = - ecl_state->get3DProperties().getDoubleGridProperty("PORV").getData(); - - Opm::GridManager gm(ecl_state->getInputGrid(), porv); - const Grid& grid = *(gm.c_grid()); - current_timestep = 0; - - // Create wells. - wells_manager.reset(new Opm::WellsManager(*ecl_state, - *schedule, - *summaryState, - current_timestep, - Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::globalCell(grid), - Opm::UgGridHelpers::cartDims(grid), - Opm::UgGridHelpers::dimensions(grid), - Opm::UgGridHelpers::cell2Faces(grid), - Opm::UgGridHelpers::beginFaceCentroids(grid), - false, - std::unordered_set() ) ); - }; - std::unique_ptr wells_manager; std::unique_ptr ecl_state; std::unique_ptr schedule; std::unique_ptr summaryState; + std::vector> well_perf_data; int current_timestep; }; @@ -132,7 +109,6 @@ BOOST_GLOBAL_FIXTURE(GlobalFixture); BOOST_AUTO_TEST_CASE(TestStandardWellInput) { const SetupTest setup_test; - const Wells* wells = setup_test.wells_manager->c_wells(); const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep); BOOST_CHECK_EQUAL( wells_ecl.size(), 2); const Opm::Well& well = wells_ecl[1]; @@ -149,36 +125,24 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) { rateConverter.reset(new RateConverterType (phaseUsage, std::vector(10, 0))); - const int pvtIdx = 0; - const int num_comp = wells->number_of_phases; + Opm::PerforationData dummy; + std::vector pdata(well.getConnections().size(), dummy); - BOOST_CHECK_THROW( StandardWell( well, -1, wells, param, *rateConverter, pvtIdx, num_comp), std::invalid_argument); - BOOST_CHECK_THROW( StandardWell( well, 4, nullptr , param, *rateConverter, pvtIdx, num_comp), std::invalid_argument); + BOOST_CHECK_THROW( StandardWell( well, -1, param, *rateConverter, 0, 3, 3, 0, 0, pdata), std::invalid_argument); } BOOST_AUTO_TEST_CASE(TestBehavoir) { const SetupTest setup_test; - const Wells* wells_struct = setup_test.wells_manager->c_wells(); const auto& wells_ecl = setup_test.schedule->getWells(setup_test.current_timestep); const int current_timestep = setup_test.current_timestep; std::vector > wells; { - const int nw = wells_struct ? (wells_struct->number_of_wells) : 0; + const int nw = wells_ecl.size(); const Opm::BlackoilModelParametersEbos param; for (int w = 0; w < nw; ++w) { - const std::string well_name(wells_struct->name[w]); - - size_t index_well = 0; - for (; index_well < wells_ecl.size(); ++index_well) { - if (well_name == wells_ecl[index_well].name()) { - break; - } - } - // we should always be able to find the well in wells_ecl - BOOST_CHECK(index_well != wells_ecl.size()); // For the conversion between the surface volume rate and resrevoir voidage rate typedef Opm::BlackOilFluidSystem FluidSystem; using RateConverterType = Opm::RateConverter:: @@ -191,11 +155,9 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { // Compute reservoir volumes for RESV controls. rateConverter.reset(new RateConverterType (phaseUsage, std::vector(10, 0))); - - const int pvtIdx = 0; - const int num_comp = wells_struct->number_of_phases; - - wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param, *rateConverter, pvtIdx, num_comp) ); + Opm::PerforationData dummy; + std::vector pdata(wells_ecl[w].getConnections().size(), dummy); + wells.emplace_back(new StandardWell(wells_ecl[w], current_timestep, param, *rateConverter, 0, 3, 3, w, 0, pdata) ); } } @@ -203,7 +165,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { { const auto& well = wells[0]; BOOST_CHECK_EQUAL(well->name(), "PROD1"); - BOOST_CHECK(well->wellType() == PRODUCER); + BOOST_CHECK(well->isProducer()); BOOST_CHECK(well->numEq == 3); BOOST_CHECK(well->numStaticWellEq== 4); } @@ -212,7 +174,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) { { const auto& well = wells[1]; BOOST_CHECK_EQUAL(well->name(), "INJE1"); - BOOST_CHECK(well->wellType() == INJECTOR); + BOOST_CHECK(well->isInjector()); BOOST_CHECK(well->numEq == 3); BOOST_CHECK(well->numStaticWellEq== 4); } diff --git a/tests/test_wellstatefullyimplicitblackoil.cpp b/tests/test_wellstatefullyimplicitblackoil.cpp index 13ea88e4d..72a440574 100644 --- a/tests/test_wellstatefullyimplicitblackoil.cpp +++ b/tests/test_wellstatefullyimplicitblackoil.cpp @@ -53,13 +53,60 @@ struct Setup , grid (es.getInputGrid()) , sched(deck, es) , st(std::chrono::system_clock::from_time_t(sched.getStartTime())) - {} + { + initWellPerfData(); + } + + void initWellPerfData() + { + const auto& wells = sched.getWells(0); + const auto& cartDims = Opm::UgGridHelpers::cartDims(*grid.c_grid()); + const int* compressed_to_cartesian = Opm::UgGridHelpers::globalCell(*grid.c_grid()); + std::vector cartesian_to_compressed(cartDims[0] * cartDims[1] * cartDims[2], -1); + for (int ii = 0; ii < Opm::UgGridHelpers::numCells(*grid.c_grid()); ++ii) { + cartesian_to_compressed[compressed_to_cartesian[ii]] = ii; + } + well_perf_data.resize(wells.size()); + int well_index = 0; + for (const auto& well : wells) { + well_perf_data[well_index].clear(); + well_perf_data[well_index].reserve(well.getConnections().size()); + for (const auto& completion : well.getConnections()) { + if (completion.state() == Opm::Connection::State::OPEN) { + const int i = completion.getI(); + const int j = completion.getJ(); + const int k = completion.getK(); + const int cart_grid_indx = i + cartDims[0] * (j + cartDims[1] * k); + const int active_index = cartesian_to_compressed[cart_grid_indx]; + if (active_index < 0) { + const std::string msg + = ("Cell with i,j,k indices " + std::to_string(i) + " " + std::to_string(j) + " " + + std::to_string(k) + " not found in grid (well = " + well.name() + ")."); + OPM_THROW(std::runtime_error, msg); + } else { + Opm::PerforationData pd; + pd.cell_index = active_index; + pd.connection_transmissibility_factor = completion.CF() * completion.wellPi(); + pd.satnum_id = completion.satTableId(); + well_perf_data[well_index].push_back(pd); + } + } else { + if (completion.state() != Opm::Connection::State::SHUT) { + OPM_THROW(std::runtime_error, + "Completion state: " << Opm::Connection::State2String(completion.state()) << " not handled"); + } + } + } + ++well_index; + } + } Opm::EclipseState es; Opm::PhaseUsage pu; Opm::GridManager grid; Opm::Schedule sched; Opm::SummaryState st; + std::vector> well_perf_data; }; namespace { @@ -72,14 +119,11 @@ namespace { std::vector(setup.grid.c_grid()->number_of_cells, 100.0*Opm::unit::barsa); - const Opm::WellsManager wmgr{setup.es, setup.sched, setup.st, timeStep, *setup.grid.c_grid()}; - - state.init(wmgr.c_wells(), cpress, setup.sched, + state.init(cpress, setup.sched, setup.sched.getWells(timeStep), - timeStep, nullptr, setup.pu); + timeStep, nullptr, setup.pu, setup.well_perf_data, setup.st); - state.initWellStateMSWell(wmgr.c_wells(), - setup.sched.getWells(timeStep), + state.initWellStateMSWell(setup.sched.getWells(timeStep), setup.pu, nullptr); return state;