From 65edfb702cd84d7ec29127488915597caddaa642 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 1 Mar 2021 15:17:52 +0100 Subject: [PATCH 1/2] Preparations for Recording Dynamic Well State Coalesce blocks with same conditions, split long lines, and apply 'const' where appropriate. While here, also tighten the "rate = 0" criterion to include denormalised numbers. --- opm/simulators/wells/BlackoilWellModel.hpp | 9 +- .../wells/BlackoilWellModel_impl.hpp | 156 ++++++++++++------ opm/simulators/wells/StandardWell_impl.hpp | 12 +- .../wells/WellStateFullyImplicitBlackoil.hpp | 100 +++++------ 4 files changed, 163 insertions(+), 114 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index f168ab290..f4e07eec8 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -268,11 +268,12 @@ namespace Opm { /// Returns true if the well was actually found and shut. bool forceShutWellByNameIfPredictionMode(const std::string& wellname, const double simulation_time); - void updateEclWell(int timeStepIdx, int well_index); - void updateEclWell(int timeStepIdx, const std::string& wname); + void updateEclWell(const int timeStepIdx, const int well_index); + void updateEclWell(const int timeStepIdx, const std::string& wname); bool hasWell(const std::string& wname); - double wellPI(int well_index) const; + double wellPI(const int well_index) const; double wellPI(const std::string& well_name) const; + protected: Simulator& ebosSimulator_; @@ -345,7 +346,7 @@ namespace Opm { WellTestState wellTestState_; std::unique_ptr guideRate_; - std::map node_pressures_; // Storing network pressures for output. + std::map node_pressures_{}; // Storing network pressures for output. // used to better efficiency of calcuation mutable BVector scaleAddRes_; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 0e2124d07..32f04ee0a 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -298,8 +298,8 @@ namespace Opm { template void BlackoilWellModel:: - beginTimeStep() { - + beginTimeStep() + { updatePerforationIntensiveQuantities(); Opm::DeferredLogger local_deferredLogger; @@ -425,17 +425,22 @@ namespace Opm { template void - BlackoilWellModel::wellTesting(const int timeStepIdx, const double simulationTime, Opm::DeferredLogger& deferred_logger) { + BlackoilWellModel::wellTesting(const int timeStepIdx, + const double simulationTime, + Opm::DeferredLogger& deferred_logger) + { const auto& wtest_config = schedule()[timeStepIdx].wtest_config(); if (wtest_config.size() != 0) { // there is a WTEST request // average B factors are required for the convergence checking of well equations // Note: this must be done on all processes, even those with // no wells needing testing, otherwise we will have locking. - std::vector< Scalar > B_avg(numComponents(), Scalar() ); + std::vector< Scalar > B_avg(numComponents(), Scalar()); computeAverageFormationFactor(B_avg); - const auto& wellsForTesting = wellTestState_.updateWells(wtest_config, wells_ecl_, simulationTime); + const auto wellsForTesting = wellTestState_ + .updateWells(wtest_config, wells_ecl_, simulationTime); + for (const auto& testWell : wellsForTesting) { const std::string& well_name = testWell.first; @@ -444,9 +449,12 @@ namespace Opm { // some preparation before the well can be used well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg); + const Well& wellEcl = schedule().getWell(well_name, timeStepIdx); double well_efficiency_factor = wellEcl.getEfficiencyFactor(); - WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx), schedule(), timeStepIdx, well_efficiency_factor); + WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx), + schedule(), timeStepIdx, well_efficiency_factor); + well->setWellEfficiencyFactor(well_efficiency_factor); well->setVFPProperties(vfp_properties_.get()); well->setGuideRate(guideRate_.get()); @@ -459,29 +467,43 @@ namespace Opm { } } + + + + // called at the end of a report step template void BlackoilWellModel:: - endReportStep() { + endReportStep() + { // Clear the communication data structures for above values. - for(auto&& pinfo : local_parallel_well_info_) + for (auto&& pinfo : local_parallel_well_info_) { pinfo->clear(); } } + + + + // called at the end of a report step template const SimulatorReportSingle& BlackoilWellModel:: lastReport() const {return last_report_; } + + + + // called at the end of a time step template void BlackoilWellModel:: - timeStepSucceeded(const double& simulationTime, const double dt) { + timeStepSucceeded(const double& simulationTime, const double dt) + { // time step is finished and we are not any more at the beginning of an report step report_step_starts_ = false; @@ -735,12 +757,13 @@ namespace Opm { if (nw > 0) { well_container.reserve(nw); + for (int w = 0; w < nw; ++w) { const Well& well_ecl = wells_ecl_[w]; const std::string& well_name = well_ecl.name(); // A new WCON keywords can re-open a well that was closed/shut due to Physical limit - if ( wellTestState_.hasWellClosed(well_name)) { + if (this->wellTestState_.hasWellClosed(well_name)) { // TODO: more checking here, to make sure this standard more specific and complete // maybe there is some WCON keywords will not open the well if (well_state_.effectiveEventsOccurred(w)) { @@ -759,9 +782,10 @@ namespace Opm { // TODO: should we do this for all kinds of closing reasons? // something like wellTestState_.hasWell(well_name)? bool wellIsStopped = false; - if ( wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::ECONOMIC) || - wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::PHYSICAL) ) { - if( well_ecl.getAutomaticShutIn() ) { + if (wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::ECONOMIC) || + wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::PHYSICAL)) + { + if (well_ecl.getAutomaticShutIn()) { // shut wells are not added to the well container well_state_.shutWell(w); continue; @@ -783,28 +807,41 @@ namespace Opm { // (prediction type) rate control is zero, then it is effectively shut. if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) { const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - auto prod_controls = well_ecl.productionControls(summaryState); + const auto prod_controls = well_ecl.productionControls(summaryState); + + auto is_zero = [](const double x) + { + return std::isfinite(x) && !std::isnormal(x); + }; + bool zero_rate_control = false; switch (prod_controls.cmode) { case Well::ProducerCMode::ORAT: - zero_rate_control = (prod_controls.oil_rate == 0.0); + zero_rate_control = is_zero(prod_controls.oil_rate); break; + case Well::ProducerCMode::WRAT: - zero_rate_control = (prod_controls.water_rate == 0.0); + zero_rate_control = is_zero(prod_controls.water_rate); break; + case Well::ProducerCMode::GRAT: - zero_rate_control = (prod_controls.gas_rate == 0.0); + zero_rate_control = is_zero(prod_controls.gas_rate); break; + case Well::ProducerCMode::LRAT: - zero_rate_control = (prod_controls.liquid_rate == 0.0); + zero_rate_control = is_zero(prod_controls.liquid_rate); break; + case Well::ProducerCMode::RESV: - zero_rate_control = (prod_controls.resv_rate == 0.0); + zero_rate_control = is_zero(prod_controls.resv_rate); break; + default: // Might still have zero rate controls, but is pressure controlled. zero_rate_control = false; + break; } + if (zero_rate_control) { // Treat as shut, do not add to container. local_deferredLogger.info(" Well shut due to zero rate control and disallowing crossflow: " + well_ecl.name()); @@ -867,10 +904,11 @@ namespace Opm { const auto& perf_data = this->well_perf_data_[wellID]; // Cater for case where local part might have no perforations. - const int pvtreg = perf_data.empty() ? - 0 : pvt_region_idx_[perf_data.front().cell_index]; + const auto pvtreg = perf_data.empty() + ? 0 : pvt_region_idx_[perf_data.front().cell_index]; + const auto& parallel_well_info = *local_parallel_well_info_[wellID]; - auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg); + const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg); return std::make_unique(this->wells_ecl_[wellID], parallel_well_info, @@ -1426,11 +1464,13 @@ namespace Opm { } } } - logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "computeWellPotentials() failed.", terminal_output_); + + logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, + "computeWellPotentials() failed.", + terminal_output_); // Store it in the well state well_state_.wellPotentials() = well_potentials; - } @@ -2536,7 +2576,7 @@ namespace Opm { template void BlackoilWellModel:: - updateEclWell(int timeStepIdx, int well_index) + updateEclWell(const int timeStepIdx, const int well_index) { const auto& schedule = this->ebosSimulator_.vanguard().schedule(); const auto& wname = this->wells_ecl_[well_index].name(); @@ -2556,23 +2596,36 @@ namespace Opm { this->prod_index_calc_[well_index].reInit(well); } + + + + template void BlackoilWellModel:: - updateEclWell(int timeStepIdx, const std::string& wname) { - auto well_iter = std::find_if( this->wells_ecl_.begin(), this->wells_ecl_.end(), [wname] (const auto& well) -> bool { return well.name() == wname;}); - if (well_iter == this->wells_ecl_.end()) - throw std::logic_error("Could not find well: " + wname); + updateEclWell(const int timeStepIdx, const std::string& wname) + { + auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(), + [&wname](const auto& well) -> bool + { + return well.name() == wname; + }); - auto well_index = std::distance( this->wells_ecl_.begin(), well_iter ); + if (well_iter == this->wells_ecl_.end()) { + throw std::logic_error { "Could not find well: " + wname }; + } + + auto well_index = std::distance(this->wells_ecl_.begin(), well_iter); this->updateEclWell(timeStepIdx, well_index); } + + template double BlackoilWellModel:: - wellPI(int well_index) const + wellPI(const int well_index) const { const auto& pu = this->phase_usage_; const auto np = this->numPhases(); @@ -2598,21 +2651,31 @@ namespace Opm { default: throw std::invalid_argument { "Unsupported preferred phase " + - std::to_string(static_cast(preferred)) - }; + std::to_string(static_cast(preferred)) + }; } } + + + + template double BlackoilWellModel:: wellPI(const std::string& well_name) const { - auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(), [&well_name](const Well& well) { return well.name() == well_name; }); - if (well_iter == this->wells_ecl_.end()) - throw std::logic_error("Could not find well: " + well_name); + auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(), + [&well_name](const Well& well) + { + return well.name() == well_name; + }); - auto well_index = std::distance( this->wells_ecl_.begin(), well_iter ); + if (well_iter == this->wells_ecl_.end()) { + throw std::logic_error { "Could not find well: " + well_name }; + } + + auto well_index = std::distance(this->wells_ecl_.begin(), well_iter); return this->wellPI(well_index); } @@ -2630,20 +2693,20 @@ namespace Opm { auto hasWellPIEvent = [this, timeStepIdx](const int well_index) -> bool { - return this->schedule()[timeStepIdx] - .wellgroup_events().hasEvent(this->wells_ecl_[well_index].name(), - ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX); + return this->schedule()[timeStepIdx].wellgroup_events() + .hasEvent(this->wells_ecl_[well_index].name(), + ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX); }; auto rescaleWellPI = [this, timeStepIdx](const int well_index, const double newWellPI) -> void { - { - const auto& wname = this->wells_ecl_[well_index].name(); - auto& schedule = this->ebosSimulator_.vanguard().schedule(); // Mutable - schedule.applyWellProdIndexScaling(wname, timeStepIdx, newWellPI); - } + const auto& wname = this->wells_ecl_[well_index].name(); + + auto& schedule = this->ebosSimulator_.vanguard().schedule(); // Mutable + schedule.applyWellProdIndexScaling(wname, timeStepIdx, newWellPI); + this->updateEclWell(timeStepIdx, well_index); }; @@ -2672,8 +2735,7 @@ namespace Opm { const auto nw = this->numLocalWells(); for (auto wellID = 0*nw; wellID < nw; ++wellID) { if (hasWellPIEvent(wellID)) { - const auto newWellPI = this->wellPI(wellID); - rescaleWellPI(wellID, newWellPI); + rescaleWellPI(wellID, this->wellPI(wellID)); } } diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 8d092924b..ed9e113d5 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -2263,9 +2263,9 @@ namespace Opm setToZero(wellPI); const auto preferred_phase = this->well_ecl_.getPreferredPhase(); - auto subsetPerfID = 0; + auto subsetPerfID = 0; - for ( const auto& perf : *this->perf_data_){ + for (const auto& perf : *this->perf_data_) { auto allPerfID = perf.ecl_index; auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double @@ -2293,13 +2293,13 @@ namespace Opm connPI += np; } - // Sum with communication in case of distributed well. + // Sum with communication in case of distributed well. const auto& comm = this->parallel_well_info_.communication(); - if (comm.size() > 1) - { + if (comm.size() > 1) { comm.sum(wellPI, np); } - assert (static_cast(subsetPerfID) == this->number_of_perforations_ && + + assert ((static_cast(subsetPerfID) == this->number_of_perforations_) && "Internal logic error in processing connections for PI/II"); } diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 8a19d89ce..eb8016a37 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -553,99 +553,85 @@ namespace Opm data::Wells res = WellState::report(pu, globalCellIdxMap); const int nw = this->numWells(); - if( nw == 0 ) return res; + if (nw == 0) { + return res; + } + const int np = pu.num_phases; - using rt = data::Rates::opt; - std::vector< rt > phs( np ); - if( pu.phase_used[Water] ) { + std::vector phs(np); + if (pu.phase_used[Water]) { phs.at( pu.phase_pos[Water] ) = rt::wat; } - if( pu.phase_used[Oil] ) { + if (pu.phase_used[Oil]) { phs.at( pu.phase_pos[Oil] ) = rt::oil; } - if( pu.phase_used[Gas] ) { + if (pu.phase_used[Gas]) { phs.at( pu.phase_pos[Gas] ) = rt::gas; } - /* 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 - * simulators, as it relies on simulator internals, details and - * representations. - */ + // This is a reference or example on **how** to convert from + // WellState to something understood by opm-common's output + // layer. It is intended to be properly implemented and + // maintained as a part of simulators, as it relies on simulator + // internals, details and representations. - for( const auto& wt : this->wellMap() ) { + for (const auto& wt : this->wellMap()) { const auto w = wt.second[ 0 ]; const auto& pwinfo = *parallel_well_info_[w]; if ((this->status_[w] == Well::Status::SHUT) || !pwinfo.isOwner()) + { continue; + } - auto& well = res.at( wt.first ); - //well.injectionControl = static_cast(this->currentInjectionControls()[ w ]); - //well.productionControl = static_cast(this->currentProductionControls()[ w ]); + auto& well = res.at(wt.first); const int well_rate_index = w * pu.num_phases; - if ( pu.phase_used[Water] ) { - well.rates.set( rt::reservoir_water, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Water]] ); + if (pu.phase_used[Water]) { + const auto i = well_rate_index + pu.phase_pos[Water]; + well.rates.set(rt::reservoir_water, this->well_reservoir_rates_[i]); + well.rates.set(rt::productivity_index_water, this->productivity_index_[i]); + well.rates.set(rt::well_potential_water, this->well_potentials_[i]); } - if ( pu.phase_used[Oil] ) { - well.rates.set( rt::reservoir_oil, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Oil]] ); + if (pu.phase_used[Oil]) { + const auto i = well_rate_index + pu.phase_pos[Oil]; + well.rates.set(rt::reservoir_oil, this->well_reservoir_rates_[i]); + well.rates.set(rt::productivity_index_oil, this->productivity_index_[i]); + well.rates.set(rt::well_potential_oil, this->well_potentials_[i]); } - if ( pu.phase_used[Gas] ) { - well.rates.set( rt::reservoir_gas, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Gas]] ); + if (pu.phase_used[Gas]) { + const auto i = well_rate_index + pu.phase_pos[Gas]; + well.rates.set(rt::reservoir_gas, this->well_reservoir_rates_[i]); + well.rates.set(rt::productivity_index_gas, this->productivity_index_[i]); + well.rates.set(rt::well_potential_gas, this->well_potentials_[i]); } - if ( pu.phase_used[Water] ) { - well.rates.set( rt::productivity_index_water, this->productivity_index_[well_rate_index + pu.phase_pos[Water]] ); + if (pu.has_solvent || pu.has_zFraction) { + well.rates.set(rt::solvent, solventWellRate(w)); } - if ( pu.phase_used[Oil] ) { - well.rates.set( rt::productivity_index_oil, this->productivity_index_[well_rate_index + pu.phase_pos[Oil]] ); + if (pu.has_polymer) { + well.rates.set(rt::polymer, polymerWellRate(w)); } - if ( pu.phase_used[Gas] ) { - well.rates.set( rt::productivity_index_gas, this->productivity_index_[well_rate_index + pu.phase_pos[Gas]] ); + if (pu.has_brine) { + well.rates.set(rt::brine, brineWellRate(w)); } - if ( pu.phase_used[Water] ) { - well.rates.set( rt::well_potential_water, this->well_potentials_[well_rate_index + pu.phase_pos[Water]] ); - } - - if ( pu.phase_used[Oil] ) { - well.rates.set( rt::well_potential_oil, this->well_potentials_[well_rate_index + pu.phase_pos[Oil]] ); - } - - if ( pu.phase_used[Gas] ) { - well.rates.set( rt::well_potential_gas, this->well_potentials_[well_rate_index + pu.phase_pos[Gas]] ); - } - - if ( pu.has_solvent || pu.has_zFraction) { - well.rates.set( rt::solvent, solventWellRate(w) ); - } - - if ( pu.has_polymer ) { - well.rates.set( rt::polymer, polymerWellRate(w) ); - } - - if ( pu.has_brine ) { - well.rates.set( rt::brine, brineWellRate(w) ); - } - - if ( is_producer_[w] ) { - well.rates.set( rt::alq, getALQ(/*wellName=*/wt.first) ); + if (is_producer_[w]) { + well.rates.set(rt::alq, getALQ(/*wellName=*/wt.first)); } else { - well.rates.set( rt::alq, 0.0 ); + well.rates.set(rt::alq, 0.0); } - well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] ); - well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] ); + well.rates.set(rt::dissolved_gas, this->well_dissolved_gas_rates_[w]); + well.rates.set(rt::vaporized_oil, this->well_vaporized_oil_rates_[w]); { auto& curr = well.current_control; From b982ad0fd2a67ec074ec5ac245f170456a0a4cb2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 1 Mar 2021 15:43:43 +0100 Subject: [PATCH 2/2] Record Dynamic Well Status in 'wellData()' Output This commit sets the 'data::Well::dynamicStatus' based on the dynamically updated 'Schedule' object (i.e., from ACTIONX and similar) and the results of well/operability testing (WECON and/or WTEST). If a well is closed due to economic limits (WECON) we still provide summary-style data at the timestep that closed the well, but omit this data at later steps until the well reopens. We add a new parameter to WellState::report() to distinguish these situations. This is in preparation of making the 'BlackoilWellModel' manage both open and shut wells alike. --- opm/simulators/wells/BlackoilWellModel.hpp | 12 ++- .../wells/BlackoilWellModel_impl.hpp | 84 ++++++++++++++++--- opm/simulators/wells/WellState.hpp | 12 ++- .../wells/WellStateFullyImplicitBlackoil.hpp | 16 ++-- tests/test_wellstatefullyimplicitblackoil.cpp | 4 +- 5 files changed, 105 insertions(+), 23 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index f4e07eec8..b41bca0e2 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -217,7 +217,12 @@ namespace Opm { Opm::data::Wells wellData() const { - auto wsrpt = well_state_.report(phase_usage_, Opm::UgGridHelpers::globalCell(grid())); + auto wsrpt = well_state_ + .report(phase_usage_, Opm::UgGridHelpers::globalCell(grid()), + [this](const int well_ndex) -> bool + { + return this->wasDynamicallyShutThisTimeStep(well_ndex); + }); this->assignWellGuideRates(wsrpt); this->assignShutConnections(wsrpt); @@ -347,6 +352,7 @@ namespace Opm { std::unique_ptr guideRate_; std::map node_pressures_{}; // Storing network pressures for output. + mutable std::unordered_set closed_this_step_{}; // used to better efficiency of calcuation mutable BVector scaleAddRes_; @@ -430,6 +436,8 @@ namespace Opm { int numPhases() const; + int reportStepIndex() const; + void assembleWellEq(const std::vector& B_avg, const double dt, Opm::DeferredLogger& deferred_logger); // some preparation work, mostly related to group control and RESV, @@ -489,6 +497,8 @@ namespace Opm { void runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger); + bool wasDynamicallyShutThisTimeStep(const int well_index) const; + void assignWellGuideRates(data::Wells& wsrpt) const; void assignShutConnections(data::Wells& wsrpt) const; void assignGroupValues(const int reportStepIdx, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 32f04ee0a..9fe6136fb 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -504,6 +504,7 @@ namespace Opm { BlackoilWellModel:: timeStepSucceeded(const double& simulationTime, const double dt) { + this->closed_this_step_.clear(); // time step is finished and we are not any more at the beginning of an report step report_step_starts_ = false; @@ -761,6 +762,20 @@ namespace Opm { for (int w = 0; w < nw; ++w) { const Well& well_ecl = wells_ecl_[w]; const std::string& well_name = well_ecl.name(); + const auto well_status = this->schedule() + .getWell(well_name, time_step).getStatus(); + + if ((well_ecl.getStatus() == Well::Status::SHUT) || + (well_status == Well::Status::SHUT)) + { + // Due to ACTIONX the well might have been closed behind our back. + if (well_ecl.getStatus() != Well::Status::SHUT) { + this->closed_this_step_.insert(well_name); + well_state_.shutWell(w); + } + + continue; + } // A new WCON keywords can re-open a well that was closed/shut due to Physical limit if (this->wellTestState_.hasWellClosed(well_name)) { @@ -796,13 +811,6 @@ namespace Opm { } } - // Due to ACTIONX the well might have been closed 'behind our back'. - const auto well_status = schedule().getWell(well_name, time_step).getStatus(); - if (well_status == Well::Status::SHUT) { - well_state_.shutWell(w); - continue; - } - // If a production well disallows crossflow and its // (prediction type) rate control is zero, then it is effectively shut. if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) { @@ -1410,8 +1418,15 @@ namespace Opm { { Opm::DeferredLogger local_deferredLogger; for (const auto& well : well_container_) { + const auto wasClosed = wellTestState.hasWellClosed(well->name()); + well->updateWellTestState(well_state_, simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger); + + if (!wasClosed && wellTestState.hasWellClosed(well->name())) { + this->closed_this_step_.insert(well->name()); + } } + Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger); if (terminal_output_) { global_deferredLogger.logMessages(); @@ -2680,6 +2695,21 @@ namespace Opm { } + + + + template + int + BlackoilWellModel:: + reportStepIndex() const + { + return std::max(this->ebosSimulator_.episodeIndex(), 0); + } + + + + + template void BlackoilWellModel:: @@ -2746,6 +2776,19 @@ namespace Opm { + template + bool + BlackoilWellModel:: + wasDynamicallyShutThisTimeStep(const int well_index) const + { + return this->closed_this_step_.find(this->wells_ecl_[well_index].name()) + != this->closed_this_step_.end(); + } + + + + + template void BlackoilWellModel:: @@ -2826,15 +2869,30 @@ namespace Opm { BlackoilWellModel:: assignShutConnections(data::Wells& wsrpt) const { + auto wellID = 0; + for (const auto& well : this->wells_ecl_) { - auto xwPos = wsrpt.find(well.name()); - if (xwPos == wsrpt.end()) { // No well results. Unexpected. - continue; + auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<> + + xwel.dynamicStatus = this->schedule() + .getWell(well.name(), this->reportStepIndex()).getStatus(); + + const auto wellIsOpen = xwel.dynamicStatus == Well::Status::OPEN; + auto skip = [wellIsOpen](const Connection& conn) + { + return wellIsOpen && (conn.state() != Connection::State::SHUT); + }; + + if (this->wellTestState_.hasWellClosed(well.name()) && + !this->wasDynamicallyShutThisTimeStep(wellID)) + { + xwel.dynamicStatus = well.getAutomaticShutIn() + ? Well::Status::SHUT : Well::Status::STOP; } - auto& xcon = xwPos->second.connections; + auto& xcon = xwel.connections; for (const auto& conn : well.getConnections()) { - if (conn.state() != Connection::State::SHUT) { + if (skip(conn)) { continue; } @@ -2845,6 +2903,8 @@ namespace Opm { xc.effective_Kh = conn.Kh(); xc.trans_factor = conn.CF(); } + + ++wellID; } } diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 0e354ea63..bd4059d26 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -207,6 +208,7 @@ namespace Opm return well_info.isOwner(); } + /// The number of wells present. int numWells() const { @@ -237,15 +239,21 @@ namespace Opm this->thp_[well_index] = 0; } - virtual data::Wells report(const PhaseUsage& pu, const int* globalCellIdxMap) const + virtual data::Wells + report(const PhaseUsage& pu, + const int* globalCellIdxMap, + const std::function& wasDynamicallyClosed) const { using rt = data::Rates::opt; data::Wells dw; for( const auto& itr : this->wellMap_ ) { const auto well_index = itr.second[ 0 ]; - if (this->status_[well_index] == Well::Status::SHUT) + if ((this->status_[well_index] == Well::Status::SHUT) && + ! wasDynamicallyClosed(well_index)) + { continue; + } const auto& pwinfo = *parallel_well_info_[well_index]; using WellT = std::remove_reference_t; diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index eb8016a37..e1f120dec 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -546,11 +547,13 @@ namespace Opm return it->second; } - - - data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override + data::Wells + report(const PhaseUsage &pu, + const int* globalCellIdxMap, + const std::function& wasDynamicallyClosed) const override { - data::Wells res = WellState::report(pu, globalCellIdxMap); + data::Wells res = + WellState::report(pu, globalCellIdxMap, wasDynamicallyClosed); const int nw = this->numWells(); if (nw == 0) { @@ -581,8 +584,9 @@ namespace Opm for (const auto& wt : this->wellMap()) { const auto w = wt.second[ 0 ]; - const auto& pwinfo = *parallel_well_info_[w]; - if ((this->status_[w] == Well::Status::SHUT) || !pwinfo.isOwner()) + if (((this->status_[w] == Well::Status::SHUT) && + ! wasDynamicallyClosed(w)) || + ! this->parallel_well_info_[w]->isOwner()) { continue; } diff --git a/tests/test_wellstatefullyimplicitblackoil.cpp b/tests/test_wellstatefullyimplicitblackoil.cpp index b480d509a..12886728f 100644 --- a/tests/test_wellstatefullyimplicitblackoil.cpp +++ b/tests/test_wellstatefullyimplicitblackoil.cpp @@ -274,7 +274,7 @@ BOOST_AUTO_TEST_CASE(Pressure) setSegPress(wells, wstate); - const auto rpt = wstate.report(setup.pu, setup.grid.c_grid()->global_cell); + const auto rpt = wstate.report(setup.pu, setup.grid.c_grid()->global_cell, [](const int){return false;}); { const auto& xw = rpt.at("INJE01"); @@ -323,7 +323,7 @@ BOOST_AUTO_TEST_CASE(Rates) setSegRates(wells, pu, wstate); - const auto rpt = wstate.report(pu, setup.grid.c_grid()->global_cell); + const auto rpt = wstate.report(pu, setup.grid.c_grid()->global_cell, [](const int){return false;}); const auto wat = pu.phase_used[Opm::BlackoilPhases::Aqua]; const auto oil = pu.phase_used[Opm::BlackoilPhases::Liquid];