diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index a21c44383..04fbc9f66 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -55,6 +55,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/utils/ParallelFileMerger.cpp opm/simulators/utils/ParallelRestart.cpp opm/simulators/wells/ALQState.cpp + opm/simulators/wells/BlackoilWellModelGeneric.cpp opm/simulators/wells/GasLiftSingleWellGeneric.cpp opm/simulators/wells/GasLiftStage2.cpp opm/simulators/wells/GlobalWellInfo.cpp diff --git a/ebos/eclwellmanager.hh b/ebos/eclwellmanager.hh index f9262f1f6..06acf48fd 100644 --- a/ebos/eclwellmanager.hh +++ b/ebos/eclwellmanager.hh @@ -568,8 +568,7 @@ public: } data::GroupAndNetworkValues - groupAndNetworkData(const int /* reportStepIdx */, - const Schedule& /* sched */) const + groupAndNetworkData(const int /* reportStepIdx */) { return {}; } diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index 883d091c9..7d370307a 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -171,7 +171,7 @@ public: const auto localWellData = simulator_.problem().wellModel().wellData(); const auto localGroupAndNetworkData = simulator_.problem().wellModel() - .groupAndNetworkData(reportStepNum, simulator_.vanguard().schedule()); + .groupAndNetworkData(reportStepNum); const auto localAquiferData = simulator_.problem().aquiferModel().aquiferData(); @@ -229,7 +229,7 @@ public: // output using eclWriter if enabled auto localWellData = simulator_.problem().wellModel().wellData(); auto localGroupAndNetworkData = simulator_.problem().wellModel() - .groupAndNetworkData(reportStepNum, simulator_.vanguard().schedule()); + .groupAndNetworkData(reportStepNum); auto localAquiferData = simulator_.problem().aquiferModel().aquiferData(); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 939853fef..4b09a5ba3 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -27,11 +27,7 @@ #include #include -#include -#include - #include -#include #include #include #include @@ -48,7 +44,6 @@ #include #include #include -#include #include #include @@ -76,6 +71,8 @@ #include +#include + namespace Opm::Properties { template @@ -90,6 +87,7 @@ namespace Opm { /// Class for handling the blackoil well model. template class BlackoilWellModel : public BaseAuxiliaryModule + , public BlackoilWellModelGeneric { public: // --------- Types --------- @@ -212,126 +210,12 @@ namespace Opm { using WellInterfacePtr = std::shared_ptr >; WellInterfacePtr well(const std::string& wellName) const; - void initFromRestartFile(const RestartValue& restartValues); - - data::GroupAndNetworkValues - groupAndNetworkData(const int reportStepIdx, const Schedule& sched) const + using BlackoilWellModelGeneric::initFromRestartFile; + void initFromRestartFile(const RestartValue& restartValues) { - auto grp_nwrk_values = ::Opm::data::GroupAndNetworkValues{}; - - this->assignGroupValues(reportStepIdx, sched, grp_nwrk_values.groupData); - this->assignNodeValues(grp_nwrk_values.nodeData); - - return grp_nwrk_values; - } - - - /* - The dynamic state of the well model is maintained with an instance - of the WellState class. Currently we have - three different wellstate instances: - - 1. The currently active wellstate is in the active_well_state_ - member. That is the state which is mutated by the simulator. - - 2. In the case timestep fails to converge and we must go back and - try again with a smaller timestep we need to recover the last - valid wellstate. This is maintained with the - last_valid_well_state_ member and the functions - commitWellState() and resetWellState(). - - 3. For the NUPCOL functionality we should either use the - currently active wellstate or a wellstate frozen at max - nupcol iterations. This is handled with the member - nupcol_well_state_ and the initNupcolWellState() function. - */ - - - /* - Immutable version of the currently active wellstate. - */ - const WellState& wellState() const - { - return this->active_wgstate_.well_state; - } - - /* - Mutable version of the currently active wellstate. - */ - WellState& wellState() - { - return this->active_wgstate_.well_state; - } - - /* - Will return the last good wellstate. This is typcially used when - initializing a new report step where the Schedule object might - have introduced new wells. The wellstate returned by - prevWellState() must have been stored with the commitWellState() - function first. - */ - const WellState& prevWellState() const - { - return this->last_valid_wgstate_.well_state; - } - - const WGState& prevWGState() const - { - return this->last_valid_wgstate_; - } - /* - Will return the currently active nupcolWellState; must initialize - the internal nupcol wellstate with initNupcolWellState() first. - */ - const WellState& nupcolWellState() const - { - return this->nupcol_wgstate_.well_state; - } - - /* - Will assign the internal member last_valid_well_state_ to the - current value of the this->active_well_state_. The state stored - with storeWellState() can then subsequently be recovered with the - resetWellState() method. - */ - void commitWGState() - { - this->last_valid_wgstate_ = this->active_wgstate_; - } - - /* - Will store a copy of the input argument well_state in the - last_valid_well_state_ member, that state can then be recovered - with a subsequent call to resetWellState(). - */ - void commitWGState(WGState wgstate) - { - this->last_valid_wgstate_ = std::move(wgstate); - } - - /* - Will update the internal variable active_well_state_ to whatever - was stored in the last_valid_well_state_ member. This function - works in pair with commitWellState() which should be called first. - */ - void resetWGState() - { - this->active_wgstate_ = this->last_valid_wgstate_; - } - - /* - Will store the current active wellstate in the nupcol_well_state_ - member. This can then be subsequently retrieved with accessor - nupcolWellState(). - */ - void updateNupcolWGState() - { - this->nupcol_wgstate_ = this->active_wgstate_; - } - - const GroupState& groupState() const - { - return this->active_wgstate_.group_state; + initFromRestartFile(restartValues, + UgGridHelpers::numCells(grid()), + param_.use_multisegment_well_); } data::Wells wellData() const @@ -343,7 +227,7 @@ namespace Opm { }); this->assignWellGuideRates(wsrpt); - this->assignShutConnections(wsrpt); + this->assignShutConnections(wsrpt, this->reportStepIndex()); return wsrpt; } @@ -365,8 +249,6 @@ namespace Opm { // Check if well equations is converged. ConvergenceReport getWellConvergence(const std::vector& B_avg, const bool checkGroupConvergence = false) const; - const PhaseUsage& phaseUsage() const { return phase_usage_; } - const SimulatorReportSingle& lastReport() const; void addWellContributions(SparseMatrixAdapter& jacobian) const @@ -386,11 +268,6 @@ namespace Opm { /// Returns true if the well was actually found and shut. bool forceShutWellByNameIfPredictionMode(const std::string& wellname, const double simulation_time); - void updateEclWells(const int timeStepIdx, const std::unordered_set& wells); - bool hasWell(const std::string& wname); - double wellPI(const int well_index) const; - double wellPI(const std::string& well_name) const; - void updatePerforationIntensiveQuantities(); // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation() // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions @@ -404,33 +281,17 @@ namespace Opm { void initPrimaryVariablesEvaluation() const; void updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls); WellInterfacePtr getWell(const std::string& well_name) const; + protected: Simulator& ebosSimulator_; - std::vector< Well > wells_ecl_{}; - std::vector< std::vector > well_perf_data_{}; - std::vector< WellProdIndexCalculator > prod_index_calc_{}; std::vector local_shut_wells_{}; - std::vector< ParallelWellInfo > parallel_well_info_; - std::vector< ParallelWellInfo* > local_parallel_well_info_; - - bool wells_active_{false}; - // a vector of all the wells. std::vector well_container_{}; - // Map from logically cartesian cell indices to compressed ones. - // Cells not in the interior are not mapped. This deactivates - // these for distributed wells and makes the distribution non-overlapping. - std::vector cartesian_to_compressed_{}; - std::vector is_cell_perforated_{}; - std::function not_on_process_{}; - - void initializeWellProdIndCalculators(); - void initializeWellPerfData(); void initializeWellState(const int timeStepIdx, const SummaryState& summaryState); @@ -452,15 +313,11 @@ namespace Opm { const ModelParameters param_; - bool terminal_output_{false}; - std::vector pvt_region_idx_{}; - PhaseUsage phase_usage_; size_t global_num_cells_{}; // the number of the cells in the local grid size_t local_num_cells_{}; double gravity_{}; std::vector depth_{}; - bool initial_step_{}; bool report_step_starts_{}; bool glift_debug = false; bool alternative_well_rate_init_{}; @@ -472,12 +329,6 @@ namespace Opm { SimulatorReportSingle last_report_{}; - WellTestState wellTestState_{}; - 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_{}; @@ -489,23 +340,10 @@ namespace Opm { const EclipseState& eclState() const { return ebosSimulator_.vanguard().eclState(); } - const Schedule& schedule() const - { return ebosSimulator_.vanguard().schedule(); } - void gliftDebug( const std::string &msg, DeferredLogger& deferred_logger) const; - /// \brief Get the wells of our partition that are not shut. - /// \param timeStepIdx The index of the time step. - /// \param[out] globalNumWells the number of wells globally. - std::vector< Well > getLocalWells(const int timeStepId) const; - - /// \brief Create the parallel well information - /// \param localWells The local wells from ECL schedule - std::vector< ParallelWellInfo* > - createLocalParallelWellInfo(const std::vector& localWells); - // compute the well fluxes and assemble them in to the reservoir equations as source terms // and in the well equations. void assemble(const int iterationIdx, @@ -521,7 +359,6 @@ namespace Opm { // xw to update Well State void recoverWellSolutionAndUpdateWellState(const BVector& x); - void updateAndCommunicateGroupData(); void updateNetworkPressures(); // setting the well_solutions_ based on well_state. @@ -549,10 +386,6 @@ namespace Opm { // The number of components in the model. int numComponents() const; - int numLocalWells() const; - - int numPhases() const; - int reportStepIndex() const; void assembleWellEq(const double dt, DeferredLogger& deferred_logger); @@ -569,91 +402,22 @@ namespace Opm { void extractLegacyDepth_(); - /// return true if wells are available in the reservoir - bool wellsActive() const; - - void setWellsActive(const bool wells_active); - - /// return true if wells are available on this process - bool localWellsActive() const; - /// upate the wellTestState related to economic limits void updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const; void wellTesting(const int timeStepIdx, const double simulationTime, DeferredLogger& deferred_logger); - // convert well data from opm-common to well state from opm-core - void loadRestartData( const data::Wells& wells, - const data::GroupAndNetworkValues& grpNwrkValues, - const PhaseUsage& phases, - const bool handle_ms_well, - WellState& state ); - - // whether there exists any multisegment well open on this process - bool anyMSWellOpenLocal() const; - - const Well& getWellEcl(const std::string& well_name) const; - - void updateGroupIndividualControls(DeferredLogger& deferred_logger, std::set& switched_groups); - void updateGroupIndividualControl(const Group& group, DeferredLogger& deferred_logger, std::set& switched_groups); - bool checkGroupConstraints(const Group& group, DeferredLogger& deferred_logger) const; - Group::ProductionCMode checkGroupProductionConstraints(const Group& group, DeferredLogger& deferred_logger) const; - Group::InjectionCMode checkGroupInjectionConstraints(const Group& group, const Phase& phase) const; - void checkGconsaleLimits(const Group& group, WellState& well_state, DeferredLogger& deferred_logger ); - - void updateGroupHigherControls(DeferredLogger& deferred_logger, std::set& switched_groups); - void checkGroupHigherConstraints(const Group& group, DeferredLogger& deferred_logger, std::set& switched_groups); - - void actionOnBrokenConstraints(const Group& group, const Group::ExceedAction& exceed_action, const Group::ProductionCMode& newControl, DeferredLogger& deferred_logger); - - void actionOnBrokenConstraints(const Group& group, const Group::InjectionCMode& newControl, const Phase& topUpPhase, DeferredLogger& deferred_logger); - - void updateWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, const WellState& wellState); - - void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent); + void setWellWsolvent(const std::string& name, double wsolvent) override; + void calcRates(const int fipnum, + const int pvtreg, + std::vector& resv_coeff) override; void runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger); - bool wasDynamicallyShutThisTimeStep(const int well_index) const; + void computeWellTemperature(); - void assignWellGuideRates(data::Wells& wsrpt) const; - void assignShutConnections(data::Wells& wsrpt) const; - void assignGroupValues(const int reportStepIdx, - const Schedule& sched, - std::map& gvalues) const; - - void assignNodeValues(std::map& gvalues) const; - - std::unordered_map - calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const; - - void assignGroupControl(const Group& group, data::GroupData& gdata) const; - data::GuideRateValue getGuideRateValues(const Well& well) const; - data::GuideRateValue getGuideRateValues(const Group& group) const; - data::GuideRateValue getGuideRateInjectionGroupValues(const Group& group) const; - void getGuideRateValues(const GuideRate::RateVector& qs, - const bool is_inj, - const std::string& wgname, - data::GuideRateValue& grval) const; - - void assignGroupGuideRates(const Group& group, - const std::unordered_map& groupGuideRates, - data::GroupData& gdata) const; - - void computeWellTemperature(); private: - GroupState& groupState() { return this->active_wgstate_.group_state; } BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& pu); - /* - The various wellState members should be accessed and modified - through the accessor functions wellState(), prevWellState(), - commitWellState(), resetWellState(), nupcolWellState() and - updateNupcolWellState(). - */ - WGState active_wgstate_; - WGState last_valid_wgstate_; - WGState nupcol_wgstate_; - }; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp new file mode 100644 index 000000000..f24e7bfd0 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -0,0 +1,1511 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace Opm { + +BlackoilWellModelGeneric:: +BlackoilWellModelGeneric(const Schedule& schedule, + const SummaryState& summaryState, + const EclipseState& eclState, + const PhaseUsage& phase_usage, + const Comm& comm) + : schedule_(schedule) + , summaryState_(summaryState) + , eclState_(eclState) + , comm_(comm) + , phase_usage_(phase_usage) + , active_wgstate_(phase_usage) + , last_valid_wgstate_(phase_usage) + , nupcol_wgstate_(phase_usage) +{ + // Create the guide rate container. + this->guideRate_ = std::make_unique(schedule_); + + const auto numProcs = comm_.size(); + this->not_on_process_ = [this, numProcs](const Well& well) { + if (numProcs == decltype(numProcs){1}) + return false; + + // Recall: false indicates NOT active! + const auto value = std::make_pair(well.name(), true); + auto candidate = std::lower_bound(this->parallel_well_info_.begin(), + this->parallel_well_info_.end(), + value); + + return (candidate == this->parallel_well_info_.end()) + || (*candidate != value); + }; +} + +int +BlackoilWellModelGeneric:: +numLocalWells() const +{ + return wells_ecl_.size(); +} + +int +BlackoilWellModelGeneric:: +numPhases() const +{ + return phase_usage_.num_phases; +} + +bool +BlackoilWellModelGeneric:: +hasWell(const std::string& wname) +{ + auto iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(), + [&wname](const Well& well) { return well.name() == wname; }); + return (iter != this->wells_ecl_.end()); +} + +bool +BlackoilWellModelGeneric:: +wellsActive() const +{ + return wells_active_; +} + +bool +BlackoilWellModelGeneric:: +localWellsActive() const +{ + return numLocalWells() > 0; +} + +bool +BlackoilWellModelGeneric:: +anyMSWellOpenLocal() const +{ + for (const auto& well : wells_ecl_) { + if (well.isMultiSegment()) { + return true; + } + } + return false; +} + +const Well& +BlackoilWellModelGeneric:: +getWellEcl(const std::string& well_name) const +{ + // finding the iterator of the well in wells_ecl + auto well_ecl = std::find_if(wells_ecl_.begin(), + wells_ecl_.end(), + [&well_name](const Well& elem)->bool { + return elem.name() == well_name; + }); + + assert(well_ecl != wells_ecl_.end()); + + return *well_ecl; +} + +void +BlackoilWellModelGeneric:: +loadRestartData(const data::Wells& rst_wells, + const data::GroupAndNetworkValues& grpNwrkValues, + const PhaseUsage& phases, + const bool handle_ms_well, + WellState& well_state) +{ + using GPMode = Group::ProductionCMode; + using GIMode = Group::InjectionCMode; + + using rt = data::Rates::opt; + const auto np = phases.num_phases; + + std::vector< rt > phs( np ); + if( phases.phase_used[BlackoilPhases::Aqua] ) { + phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat; + } + + if( phases.phase_used[BlackoilPhases::Liquid] ) { + phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil; + } + + if( phases.phase_used[BlackoilPhases::Vapour] ) { + phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; + } + + for( const auto& wm : well_state.wellMap() ) { + const auto well_index = wm.second[ 0 ]; + const auto& rst_well = rst_wells.at( wm.first ); + well_state.update_bhp(well_index, rst_well.bhp); + well_state.update_temperature(well_index, rst_well.temperature); + + if (rst_well.current_control.isProducer) { + well_state.currentProductionControl(well_index, rst_well.current_control.prod); + } + else { + well_state.currentInjectionControl(well_index, rst_well.current_control.inj); + } + + for( size_t i = 0; i < phs.size(); ++i ) { + assert( rst_well.rates.has( phs[ i ] ) ); + well_state.wellRates(well_index)[i] = rst_well.rates.get(phs[i]); + } + + auto& perf_pressure = well_state.perfPress(well_index); + auto& perf_rates = well_state.perfRates(well_index); + auto& perf_phase_rates = well_state.perfPhaseRates(well_index); + const auto& perf_data = this->well_perf_data_[well_index]; + + for (std::size_t perf_index = 0; perf_index < perf_data.size(); perf_index++) { + const auto& pd = perf_data[perf_index]; + const auto& rst_connection = rst_well.connections[pd.ecl_index]; + perf_pressure[perf_index] = rst_connection.pressure; + perf_rates[perf_index] = rst_connection.reservoir_rate; + for (int phase_index = 0; phase_index < np; ++phase_index) + perf_phase_rates[perf_index*np + phase_index] = rst_connection.rates.get(phs[phase_index]); + } + + if (handle_ms_well && !rst_well.segments.empty()) { + // we need the well_ecl_ information + const std::string& well_name = wm.first; + const Well& well_ecl = getWellEcl(well_name); + + const WellSegments& segment_set = well_ecl.getSegments(); + + const auto& rst_segments = rst_well.segments; + + // \Note: eventually we need to handle the situations that some segments are shut + assert(0u + segment_set.size() == rst_segments.size()); + + auto& segments = well_state.segments(well_index); + auto& segment_pressure = segments.pressure; + auto& segment_rates = segments.rates; + for (const auto& rst_segment : rst_segments) { + const int segment_index = segment_set.segmentNumberToIndex(rst_segment.first); + + // recovering segment rates and pressure from the restart values + const auto pres_idx = data::SegmentPressures::Value::Pressure; + segment_pressure[segment_index] = rst_segment.second.pressures[pres_idx]; + + const auto& rst_segment_rates = rst_segment.second.rates; + for (int p = 0; p < np; ++p) { + segment_rates[segment_index * np + p] = rst_segment_rates.get(phs[p]); + } + } + } + } + + for (const auto& [group, value] : grpNwrkValues.groupData) { + const auto cpc = value.currentControl.currentProdConstraint; + const auto cgi = value.currentControl.currentGasInjectionConstraint; + const auto cwi = value.currentControl.currentWaterInjectionConstraint; + + if (cpc != GPMode::NONE) { + this->groupState().production_control(group, cpc); + } + + if (cgi != GIMode::NONE) { + this->groupState().injection_control(group, Phase::GAS, cgi); + } + + if (cwi != GIMode::NONE) { + this->groupState().injection_control(group, Phase::WATER, cwi); + } + } +} + +void +BlackoilWellModelGeneric:: +initFromRestartFile(const RestartValue& restartValues, + const size_t numCells, + bool handle_ms_well) +{ + // The restart step value is used to identify wells present at the given + // time step. Wells that are added at the same time step as RESTART is initiated + // will not be present in a restart file. Use the previous time step to retrieve + // wells that have information written to the restart file. + const int report_step = std::max(eclState_.getInitConfig().getRestartStep() - 1, 0); + // wells_ecl_ should only contain wells on this processor. + wells_ecl_ = getLocalWells(report_step); + local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_); + + this->initializeWellProdIndCalculators(); + initializeWellPerfData(); + + const int nw = wells_ecl_.size(); + if (nw > 0) { + handle_ms_well &= anyMSWellOpenLocal(); + this->wellState().resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, well_perf_data_, summaryState_); // Resize for restart step + loadRestartData(restartValues.wells, restartValues.grp_nwrk, phase_usage_, handle_ms_well, this->wellState()); + } + + this->commitWGState(); + initial_step_ = false; +} + +void +BlackoilWellModelGeneric:: +setWellsActive(const bool wells_active) +{ + wells_active_ = wells_active; +} + +std::vector +BlackoilWellModelGeneric:: +getLocalWells(const int timeStepIdx) const +{ + auto w = schedule().getWells(timeStepIdx); + w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end()); + return w; +} + +std::vector +BlackoilWellModelGeneric:: +createLocalParallelWellInfo(const std::vector& wells) +{ + std::vector local_parallel_well_info; + local_parallel_well_info.reserve(wells.size()); + for (const auto& well : wells) + { + auto wellPair = std::make_pair(well.name(), true); + auto pwell = std::lower_bound(parallel_well_info_.begin(), + parallel_well_info_.end(), + wellPair); + assert(pwell != parallel_well_info_.end() && + *pwell == wellPair); + local_parallel_well_info.push_back(&(*pwell)); + } + return local_parallel_well_info; +} + +void +BlackoilWellModelGeneric:: +initializeWellProdIndCalculators() +{ + this->prod_index_calc_.clear(); + this->prod_index_calc_.reserve(this->wells_ecl_.size()); + for (const auto& well : this->wells_ecl_) { + this->prod_index_calc_.emplace_back(well); + } +} + +void +BlackoilWellModelGeneric:: +initializeWellPerfData() +{ + well_perf_data_.resize(wells_ecl_.size()); + int well_index = 0; + for (const auto& well : wells_ecl_) { + int completion_index = 0; + // INVALID_ECL_INDEX marks no above perf available + int completion_index_above = ParallelWellInfo::INVALID_ECL_INDEX; + well_perf_data_[well_index].clear(); + well_perf_data_[well_index].reserve(well.getConnections().size()); + CheckDistributedWellConnections checker(well, *local_parallel_well_info_[well_index]); + bool hasFirstPerforation = false; + bool firstOpenCompletion = true; + auto& parallelWellInfo = *local_parallel_well_info_[well_index]; + parallelWellInfo.beginReset(); + + for (const auto& completion : well.getConnections()) { + const int active_index = + cartesian_to_compressed_[completion.global_index()]; + if (completion.state() == Connection::State::OPEN) { + if (active_index >= 0) { + if (firstOpenCompletion) + { + hasFirstPerforation = true; + } + checker.connectionFound(completion_index); + PerforationData pd; + pd.cell_index = active_index; + pd.connection_transmissibility_factor = completion.CF(); + pd.satnum_id = completion.satTableId(); + pd.ecl_index = completion_index; + well_perf_data_[well_index].push_back(pd); + parallelWellInfo.pushBackEclIndex(completion_index_above, + completion_index); + } + firstOpenCompletion = false; + // Next time this index is the one above as each open completion is + // is stored somehwere. + completion_index_above = completion_index; + } else { + checker.connectionFound(completion_index); + if (completion.state() != Connection::State::SHUT) { + OPM_THROW(std::runtime_error, + "Completion state: " << Connection::State2String(completion.state()) << " not handled"); + } + } + // Note: we rely on the connections being filtered! I.e. there are only connections + // to active cells in the global grid. + ++completion_index; + } + parallelWellInfo.endReset(); + checker.checkAllConnectionsFound(); + parallelWellInfo.communicateFirstPerforation(hasFirstPerforation); + ++well_index; + } +} + +bool +BlackoilWellModelGeneric:: +checkGroupConstraints(const Group& group, + const int reportStepIdx, + DeferredLogger& deferred_logger) const +{ + if (group.isInjectionGroup()) { + const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; + for (Phase phase : all) { + if (!group.hasInjectionControl(phase)) { + continue; + } + Group::InjectionCMode newControl = checkGroupInjectionConstraints(group, reportStepIdx, phase); + if (newControl != Group::InjectionCMode::NONE) { + return true; + } + } + } + if (group.isProductionGroup()) { + Group::ProductionCMode newControl = checkGroupProductionConstraints(group, reportStepIdx, deferred_logger); + if (newControl != Group::ProductionCMode::NONE) + { + return true; + } + } + + // call recursively down the group hiearchy + bool violated = false; + for (const std::string& groupName : group.groups()) { + violated = violated || checkGroupConstraints( schedule().getGroup(groupName, reportStepIdx), reportStepIdx, deferred_logger); + } + return violated; +} + +Group::InjectionCMode +BlackoilWellModelGeneric:: +checkGroupInjectionConstraints(const Group& group, + const int reportStepIdx, + const Phase& phase) const +{ + const auto& well_state = this->wellState(); + + int phasePos; + if (phase == Phase::GAS && phase_usage_.phase_used[BlackoilPhases::Vapour] ) + phasePos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + else if (phase == Phase::OIL && phase_usage_.phase_used[BlackoilPhases::Liquid]) + phasePos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; + else if (phase == Phase::WATER && phase_usage_.phase_used[BlackoilPhases::Aqua] ) + phasePos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + else + OPM_THROW(std::runtime_error, "Unknown phase" ); + + const auto& controls = group.injectionControls(phase, summaryState_); + auto currentControl = this->groupState().injection_control(group.name(), phase); + + if (controls.has_control(Group::InjectionCMode::RATE)) + { + if (currentControl != Group::InjectionCMode::RATE) + { + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.surface_max_rate < current_rate) { + return Group::InjectionCMode::RATE; + } + } + } + if (controls.has_control(Group::InjectionCMode::RESV)) + { + if (currentControl != Group::InjectionCMode::RESV) + { + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true); + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.resv_max_rate < current_rate) { + return Group::InjectionCMode::RESV; + } + } + } + if (controls.has_control(Group::InjectionCMode::REIN)) + { + if (currentControl != Group::InjectionCMode::REIN) + { + double production_Rate = 0.0; + const Group& groupRein = schedule().getGroup(controls.reinj_group, reportStepIdx); + production_Rate += WellGroupHelpers::sumWellRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false); + + // sum over all nodes + production_Rate = comm_.sum(production_Rate); + + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.target_reinj_fraction*production_Rate < current_rate) { + return Group::InjectionCMode::REIN; + } + } + } + if (controls.has_control(Group::InjectionCMode::VREP)) + { + if (currentControl != Group::InjectionCMode::VREP) + { + double voidage_rate = 0.0; + const Group& groupVoidage = schedule().getGroup(controls.voidage_group, reportStepIdx); + voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false); + voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false); + voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false); + + // sum over all nodes + voidage_rate = comm_.sum(voidage_rate); + + double total_rate = 0.0; + total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true); + total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true); + total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true); + + // sum over all nodes + total_rate = comm_.sum(total_rate); + + if (controls.target_void_fraction*voidage_rate < total_rate) { + return Group::InjectionCMode::VREP; + } + } + } + return Group::InjectionCMode::NONE; +} + +Group::ProductionCMode +BlackoilWellModelGeneric:: +checkGroupProductionConstraints(const Group& group, + const int reportStepIdx, + DeferredLogger& deferred_logger) const +{ + const auto& well_state = this->wellState(); + + const auto controls = group.productionControls(summaryState_); + const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name()); + + if (group.has_control(Group::ProductionCMode::ORAT)) + { + if (currentControl != Group::ProductionCMode::ORAT) + { + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.oil_target < current_rate ) { + return Group::ProductionCMode::ORAT; + } + } + } + + if (group.has_control(Group::ProductionCMode::WRAT)) + { + if (currentControl != Group::ProductionCMode::WRAT) + { + + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.water_target < current_rate ) { + return Group::ProductionCMode::WRAT; + } + } + } + if (group.has_control(Group::ProductionCMode::GRAT)) + { + if (currentControl != Group::ProductionCMode::GRAT) + { + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + if (controls.gas_target < current_rate ) { + return Group::ProductionCMode::GRAT; + } + } + } + if (group.has_control(Group::ProductionCMode::LRAT)) + { + if (currentControl != Group::ProductionCMode::LRAT) + { + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false); + current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.liquid_target < current_rate ) { + return Group::ProductionCMode::LRAT; + } + } + } + + if (group.has_control(Group::ProductionCMode::CRAT)) + { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "CRAT control for production groups not implemented" , deferred_logger); + } + if (group.has_control(Group::ProductionCMode::RESV)) + { + if (currentControl != Group::ProductionCMode::RESV) + { + double current_rate = 0.0; + current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true); + current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true); + current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true); + + // sum over all nodes + current_rate = comm_.sum(current_rate); + + if (controls.resv_target < current_rate ) { + return Group::ProductionCMode::RESV; + } + } + } + if (group.has_control(Group::ProductionCMode::PRBL)) + { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "PRBL control for production groups not implemented", deferred_logger); + } + return Group::ProductionCMode::NONE; +} + +void +BlackoilWellModelGeneric:: +checkGconsaleLimits(const Group& group, + WellState& well_state, + const int reportStepIdx, + DeferredLogger& deferred_logger) +{ + // call recursively down the group hiearchy + for (const std::string& groupName : group.groups()) { + checkGconsaleLimits( schedule().getGroup(groupName, reportStepIdx), well_state, reportStepIdx, deferred_logger); + } + + // only for groups with gas injection controls + if (!group.hasInjectionControl(Phase::GAS)) { + return; + } + + // check if gconsale is used for this group + if (!schedule()[reportStepIdx].gconsale().has(group.name())) + return; + + std::ostringstream ss; + + const auto& gconsale = schedule()[reportStepIdx].gconsale().get(group.name(), summaryState_); + const Group::ProductionCMode& oldProductionControl = this->groupState().production_control(group.name()); + + int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + double production_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/false); + double injection_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/true); + + // sum over all nodes + injection_rate = comm_.sum(injection_rate); + production_rate = comm_.sum(production_rate); + + double sales_rate = production_rate - injection_rate; + double production_target = gconsale.sales_target + injection_rate; + + // add import rate and substract consumption rate for group for gas + if (schedule()[reportStepIdx].gconsump().has(group.name())) { + const auto& gconsump = schedule()[reportStepIdx].gconsump().get(group.name(), summaryState_); + if (phase_usage_.phase_used[BlackoilPhases::Vapour]) { + sales_rate += gconsump.import_rate; + sales_rate -= gconsump.consumption_rate; + production_target -= gconsump.import_rate; + production_target += gconsump.consumption_rate; + } + } + + if (sales_rate > gconsale.max_sales_rate) { + switch(gconsale.max_proc) { + case GConSale::MaxProcedure::NONE: { + if (oldProductionControl != Group::ProductionCMode::GRAT && oldProductionControl != Group::ProductionCMode::NONE) { + ss << "Group sales exceed maximum limit, but the action is NONE for " + group.name() + ". Nothing happens"; + } + break; + } + case GConSale::MaxProcedure::CON: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit CON not implemented", deferred_logger); + break; + } + case GConSale::MaxProcedure::CON_P: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit CON_P not implemented", deferred_logger); + break; + } + case GConSale::MaxProcedure::WELL: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit WELL not implemented", deferred_logger); + break; + } + case GConSale::MaxProcedure::PLUG: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit PLUG not implemented", deferred_logger); + break; + } + case GConSale::MaxProcedure::MAXR: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit MAXR not implemented", deferred_logger); + break; + } + case GConSale::MaxProcedure::END: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit END not implemented", deferred_logger); + break; + } + case GConSale::MaxProcedure::RATE: { + this->groupState().production_control(group.name(), Group::ProductionCMode::GRAT); + ss << "Maximum GCONSALE limit violated for " << group.name() << ". The group is switched from "; + ss << Group::ProductionCMode2String(oldProductionControl) << " to " << Group::ProductionCMode2String(Group::ProductionCMode::GRAT); + ss << " and limited by the maximum sales rate after consumption and import are considered" ; + this->groupState().update_grat_sales_target(group.name(), production_target); + break; + } + default: + throw("Invalid procedure for maximum rate limit selected for group" + group.name()); + } + } + if (sales_rate < gconsale.min_sales_rate) { + const Group::ProductionCMode& currentProductionControl = this->groupState().production_control(group.name()); + if ( currentProductionControl == Group::ProductionCMode::GRAT ) { + ss << "Group " + group.name() + " has sale rate less then minimum permitted value and is under GRAT control. \n"; + ss << "The GRAT is increased to meet the sales minimum rate. \n"; + this->groupState().update_grat_sales_target(group.name(), production_target); + //} else if () {//TODO add action for WGASPROD + //} else if () {//TODO add action for drilling queue + } else { + ss << "Group " + group.name() + " has sale rate less then minimum permitted value but cannot increase the group production rate \n"; + ss << "or adjust gas production using WGASPROD or drill new wells to meet the sales target. \n"; + ss << "Note that WGASPROD and drilling queues are not implemented in Flow. No action is taken. \n "; + } + } + if (gconsale.sales_target < 0.0) { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + " has sale rate target less then zero. Not implemented in Flow" , deferred_logger); + } + + if (!ss.str().empty() && comm_.rank() == 0) + deferred_logger.info(ss.str()); +} + +void +BlackoilWellModelGeneric:: +checkGroupHigherConstraints(const Group& group, + DeferredLogger& deferred_logger, + const int reportStepIdx, + std::set& switched_groups) +{ + // Set up coefficients for RESV <-> surface rate conversion. + // Use the pvtRegionIdx from the top cell of the first well. + // TODO fix this! + // This is only used for converting RESV rates. + // What is the proper approach? + const int fipnum = 0; + int pvtreg = well_perf_data_.empty() || well_perf_data_[0].empty() + ? pvt_region_idx_[0] + : pvt_region_idx_[well_perf_data_[0][0].cell_index]; + + if ( comm_.size() > 1) + { + // Just like in the sequential case the pvtregion is determined + // by the first cell of the first well. What is the first well + // is decided by the order in the Schedule using Well::seqIndex() + int firstWellIndex = well_perf_data_.empty() ? + std::numeric_limits::max() : wells_ecl_[0].seqIndex(); + auto regIndexPair = std::make_pair(pvtreg, firstWellIndex); + std::vector pairs(comm_.size()); + comm_.allgather(®IndexPair, 1, pairs.data()); + pvtreg = std::min_element(pairs.begin(), pairs.end(), + [](const auto& p1, const auto& p2){ return p1.second < p2.second;}) + ->first; + } + + std::vector resv_coeff(phase_usage_.num_phases, 0.0); + calcRates(fipnum, pvtreg, resv_coeff); + + std::vector rates(phase_usage_.num_phases, 0.0); + + const bool skip = switched_groups.count(group.name()) || group.name() == "FIELD"; + + if (!skip && group.isInjectionGroup()) { + // Obtain rates for group. + for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) { + const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ true); + // Sum over all processes + rates[phasePos] = comm_.sum(local_current_rate); + } + const Phase all[] = { Phase::WATER, Phase::OIL, Phase::GAS }; + for (Phase phase : all) { + // Check higher up only if under individual (not FLD) control. + auto currentControl = this->groupState().injection_control(group.name(), phase); + if (currentControl != Group::InjectionCMode::FLD && group.injectionGroupControlAvailable(phase)) { + const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx); + const std::pair changed = WellGroupHelpers::checkGroupConstraintsInj( + group.name(), + group.parent(), + parentGroup, + this->wellState(), + this->groupState(), + reportStepIdx, + guideRate_.get(), + rates.data(), + phase, + phase_usage_, + group.getGroupEfficiencyFactor(), + schedule(), + summaryState_, + resv_coeff, + deferred_logger); + if (changed.first) { + switched_groups.insert(group.name()); + actionOnBrokenConstraints(group, Group::InjectionCMode::FLD, phase, deferred_logger); + } + } + } + } + + if (!skip && group.isProductionGroup()) { + // Obtain rates for group. + for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) { + const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ false); + // Sum over all processes + rates[phasePos] = -comm_.sum(local_current_rate); + } + // Check higher up only if under individual (not FLD) control. + const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name()); + if (currentControl != Group::ProductionCMode::FLD && group.productionGroupControlAvailable()) { + const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx); + const std::pair changed = WellGroupHelpers::checkGroupConstraintsProd( + group.name(), + group.parent(), + parentGroup, + this->wellState(), + this->groupState(), + reportStepIdx, + guideRate_.get(), + rates.data(), + phase_usage_, + group.getGroupEfficiencyFactor(), + schedule(), + summaryState_, + resv_coeff, + deferred_logger); + if (changed.first) { + switched_groups.insert(group.name()); + const auto exceed_action = group.productionControls(summaryState_).exceed_action; + actionOnBrokenConstraints(group, exceed_action, Group::ProductionCMode::FLD, deferred_logger); + } + } + } + + // call recursively down the group hiearchy + for (const std::string& groupName : group.groups()) { + checkGroupHigherConstraints( schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx, switched_groups); + } +} + +void +BlackoilWellModelGeneric:: +updateGroupIndividualControl(const Group& group, + DeferredLogger& deferred_logger, + const int reportStepIdx, + std::set& switched_groups) +{ + const bool skip = switched_groups.count(group.name()); + if (!skip && group.isInjectionGroup()) + { + const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; + for (Phase phase : all) { + if (!group.hasInjectionControl(phase)) { + continue; + } + Group::InjectionCMode newControl = checkGroupInjectionConstraints(group, reportStepIdx, phase); + if (newControl != Group::InjectionCMode::NONE) + { + switched_groups.insert(group.name()); + actionOnBrokenConstraints(group, newControl, phase, deferred_logger); + } + } + } + if (!skip && group.isProductionGroup()) { + Group::ProductionCMode newControl = checkGroupProductionConstraints(group, reportStepIdx, deferred_logger); + const auto controls = group.productionControls(summaryState_); + if (newControl != Group::ProductionCMode::NONE) + { + switched_groups.insert(group.name()); + actionOnBrokenConstraints(group, controls.exceed_action, newControl, deferred_logger); + } + } + + // call recursively down the group hiearchy + for (const std::string& groupName : group.groups()) { + updateGroupIndividualControl( schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx, switched_groups); + } +} + +void +BlackoilWellModelGeneric:: +updateGroupIndividualControls(DeferredLogger& deferred_logger, + std::set& switched_groups, + const int reportStepIdx, + const int iterationIdx) +{ + const int nupcol = schedule()[reportStepIdx].nupcol(); + // don't switch group control when iterationIdx > nupcol + // to avoid oscilations between group controls + if (iterationIdx > nupcol) + return; + + const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); + updateGroupIndividualControl(fieldGroup, deferred_logger, + reportStepIdx, switched_groups); +} + +void +BlackoilWellModelGeneric:: +updateGroupHigherControls(DeferredLogger& deferred_logger, + const int reportStepIdx, + std::set& switched_groups) +{ + const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); + checkGroupHigherConstraints(fieldGroup, deferred_logger, reportStepIdx, switched_groups); +} + +void +BlackoilWellModelGeneric:: +actionOnBrokenConstraints(const Group& group, + const Group::ExceedAction& exceed_action, + const Group::ProductionCMode& newControl, + DeferredLogger& deferred_logger) +{ + const Group::ProductionCMode oldControl = this->groupState().production_control(group.name()); + + std::ostringstream ss; + + switch(exceed_action) { + case Group::ExceedAction::NONE: { + if (oldControl != newControl && oldControl != Group::ProductionCMode::NONE) { + ss << "Group production exceed action is NONE for group " + group.name() + ". Nothing happens."; + } + break; + } + case Group::ExceedAction::CON: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit CON not implemented", deferred_logger); + break; + } + case Group::ExceedAction::CON_PLUS: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit CON_PLUS not implemented", deferred_logger); + break; + } + case Group::ExceedAction::WELL: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit WELL not implemented", deferred_logger); + break; + } + case Group::ExceedAction::PLUG: { + OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit PLUG not implemented", deferred_logger); + break; + } + case Group::ExceedAction::RATE: { + if (oldControl != newControl) { + this->groupState().production_control(group.name(), newControl); + ss << "Switching production control mode for group "<< group.name() + << " from " << Group::ProductionCMode2String(oldControl) + << " to " << Group::ProductionCMode2String(newControl); + } + break; + } + default: + throw("Invalid procedure for maximum rate limit selected for group" + group.name()); + } + + auto cc = Dune::MPIHelper::getCollectiveCommunication(); + if (!ss.str().empty() && cc.rank() == 0) + deferred_logger.info(ss.str()); +} + +void +BlackoilWellModelGeneric:: +actionOnBrokenConstraints(const Group& group, + const Group::InjectionCMode& newControl, + const Phase& controlPhase, + DeferredLogger& deferred_logger) +{ + auto oldControl = this->groupState().injection_control(group.name(), controlPhase); + + std::ostringstream ss; + if (oldControl != newControl) { + const std::string from = Group::InjectionCMode2String(oldControl); + ss << "Switching injection control mode for group "<< group.name() + << " from " << Group::InjectionCMode2String(oldControl) + << " to " << Group::InjectionCMode2String(newControl); + this->groupState().injection_control(group.name(), controlPhase, newControl); + } + auto cc = Dune::MPIHelper::getCollectiveCommunication(); + if (!ss.str().empty() && cc.rank() == 0) + deferred_logger.info(ss.str()); +} + +void +BlackoilWellModelGeneric:: +updateEclWells(const int timeStepIdx, + const std::unordered_set& wells) +{ + for (const auto& wname : wells) { + 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()) { + auto well_index = std::distance( this->wells_ecl_.begin(), well_iter ); + this->wells_ecl_[well_index] = schedule_.getWell(wname, timeStepIdx); + + const auto& well = this->wells_ecl_[well_index]; + auto& pd = this->well_perf_data_[well_index]; + auto pdIter = pd.begin(); + for (const auto& conn : well.getConnections()) { + if (conn.state() != Connection::State::SHUT) { + pdIter->connection_transmissibility_factor = conn.CF(); + ++pdIter; + } + } + this->wellState().updateStatus(well_index, well.getStatus()); + this->wellState().resetConnectionTransFactors(well_index, pd); + this->prod_index_calc_[well_index].reInit(well); + } + } +} + +double +BlackoilWellModelGeneric:: +wellPI(const int well_index) const +{ + const auto& pu = this->phase_usage_; + const auto& pi = this->wellState().productivityIndex(well_index); + + const auto preferred = this->wells_ecl_[well_index].getPreferredPhase(); + switch (preferred) { // Should really have LIQUID = OIL + WATER here too... + case Phase::WATER: + return pu.phase_used[BlackoilPhases::PhaseIndex::Aqua] + ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Aqua]] + : 0.0; + + case Phase::OIL: + return pu.phase_used[BlackoilPhases::PhaseIndex::Liquid] + ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Liquid]] + : 0.0; + + case Phase::GAS: + return pu.phase_used[BlackoilPhases::PhaseIndex::Vapour] + ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Vapour]] + : 0.0; + + default: + throw std::invalid_argument { + "Unsupported preferred phase " + + std::to_string(static_cast(preferred)) + }; + } +} + +double +BlackoilWellModelGeneric:: +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_index = std::distance(this->wells_ecl_.begin(), well_iter); + return this->wellPI(well_index); +} + +bool +BlackoilWellModelGeneric:: +wasDynamicallyShutThisTimeStep(const int well_index) const +{ + return this->closed_this_step_.find(this->wells_ecl_[well_index].name()) != + this->closed_this_step_.end(); +} + +void +BlackoilWellModelGeneric:: +updateWsolvent(const Group& group, + const int reportStepIdx, + const WellState& wellState) +{ + for (const std::string& groupName : group.groups()) { + const Group& groupTmp = schedule_.getGroup(groupName, reportStepIdx); + updateWsolvent(groupTmp, reportStepIdx, wellState); + } + + if (group.isProductionGroup()) + return; + + auto currentGroupControl = this->groupState().injection_control(group.name(), Phase::GAS); + if( currentGroupControl == Group::InjectionCMode::REIN ) { + int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const auto& controls = group.injectionControls(Phase::GAS, summaryState_); + const Group& groupRein = schedule_.getGroup(controls.reinj_group, reportStepIdx); + double gasProductionRate = WellGroupHelpers::sumWellRates(groupRein, schedule_, wellState, reportStepIdx, gasPos, /*isInjector*/false); + double solventProductionRate = WellGroupHelpers::sumSolventRates(groupRein, schedule_, wellState, reportStepIdx, /*isInjector*/false); + + solventProductionRate = comm_.sum(solventProductionRate); + gasProductionRate = comm_.sum(gasProductionRate); + + double wsolvent = 0.0; + if (std::abs(gasProductionRate) > 1e-6) + wsolvent = solventProductionRate / gasProductionRate; + + setWsolvent(group, reportStepIdx, wsolvent); + } +} + +void +BlackoilWellModelGeneric:: +setWsolvent(const Group& group, + const int reportStepIdx, + double wsolvent) +{ + for (const std::string& groupName : group.groups()) { + const Group& groupTmp = schedule_.getGroup(groupName, reportStepIdx); + setWsolvent(groupTmp, reportStepIdx, wsolvent); + } + + for (const std::string& wellName : group.wells()) { + const auto& wellTmp = schedule_.getWell(wellName, reportStepIdx); + if (wellTmp.getStatus() == Well::Status::SHUT) + continue; + + this->setWellWsolvent(wellName, wsolvent); + } +} + +data::GuideRateValue +BlackoilWellModelGeneric:: +getGuideRateValues(const Well& well) const +{ + auto grval = data::GuideRateValue{}; + + assert (this->guideRate_ != nullptr); + + const auto& wname = well.name(); + if (!this->wellState().hasWellRates(wname)) { + // No flow rates for 'wname' -- might be before well comes + // online (e.g., for the initial condition before simulation + // starts). + return grval; + } + + if (!this->guideRate_->has(wname)) { + // No guiderates exist for 'wname'. + return grval; + } + + const auto qs = WellGroupHelpers:: + getWellRateVector(this->wellState(), this->phase_usage_, wname); + + this->getGuideRateValues(qs, well.isInjector(), wname, grval); + + return grval; +} + +void +BlackoilWellModelGeneric:: +getGuideRateValues(const GuideRate::RateVector& qs, + const bool is_inj, + const std::string& wgname, + data::GuideRateValue& grval) const +{ + auto getGR = [this, &wgname, &qs](const GuideRateModel::Target t) + { + return this->guideRate_->get(wgname, t, qs); + }; + + // Note: GuideRate does currently (2020-07-20) not support Target::RES. + grval.set(data::GuideRateValue::Item::Gas, + getGR(GuideRateModel::Target::GAS)); + + grval.set(data::GuideRateValue::Item::Water, + getGR(GuideRateModel::Target::WAT)); + + if (!is_inj) { + // Producer. Extract "all" guiderate values. + grval.set(data::GuideRateValue::Item::Oil, + getGR(GuideRateModel::Target::OIL)); + } +} + +data::GuideRateValue +BlackoilWellModelGeneric:: +getGuideRateValues(const Group& group) const +{ + auto grval = data::GuideRateValue{}; + + assert (this->guideRate_ != nullptr); + + const auto& gname = group.name(); + + if (!this->groupState().has_production_rates(gname)) { + // No flow rates for production group 'gname' -- might be before group comes + // online (e.g., for the initial condition before simulation + // starts). + return grval; + } + + if (!this->guideRate_->has(gname)) { + // No guiderates exist for 'gname'. + return grval; + } + + const auto qs = WellGroupHelpers::getProductionGroupRateVector(this->groupState(), this->phase_usage_, gname); + + const auto is_inj = false; // This procedure only applies to G*PGR. + this->getGuideRateValues(qs, is_inj, gname, grval); + + return grval; +} + +data::GuideRateValue +BlackoilWellModelGeneric:: +getGuideRateInjectionGroupValues(const Group& group) const +{ + auto grval = data::GuideRateValue{}; + + assert (this->guideRate_ != nullptr); + + const auto& gname = group.name(); + if (this->guideRate_->has(gname, Phase::GAS)) { + grval.set(data::GuideRateValue::Item::Gas, + this->guideRate_->get(gname, Phase::GAS)); + } + if (this->guideRate_->has(gname, Phase::WATER)) { + grval.set(data::GuideRateValue::Item::Water, + this->guideRate_->get(gname, Phase::WATER)); + } + return grval; +} + +void +BlackoilWellModelGeneric:: +assignWellGuideRates(data::Wells& wsrpt) const +{ + for (const auto& well : this->wells_ecl_) { + auto xwPos = wsrpt.find(well.name()); + if (xwPos == wsrpt.end()) { // No well results. Unexpected. + continue; + } + + xwPos->second.guide_rates = this->getGuideRateValues(well); + } +} + +void +BlackoilWellModelGeneric:: +assignShutConnections(data::Wells& wsrpt, + const int reportStepIndex) const +{ + auto wellID = 0; + + for (const auto& well : this->wells_ecl_) { + auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<> + + xwel.dynamicStatus = this->schedule() + .getWell(well.name(), 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 = xwel.connections; + for (const auto& conn : well.getConnections()) { + if (skip(conn)) { + continue; + } + + auto& xc = xcon.emplace_back(); + xc.index = conn.global_index(); + xc.pressure = xc.reservoir_rate = 0.0; + + xc.effective_Kh = conn.Kh(); + xc.trans_factor = conn.CF(); + } + + ++wellID; + } +} + +std::unordered_map +BlackoilWellModelGeneric:: +calculateAllGroupGuiderates(const int reportStepIdx) const +{ + auto gr = std::unordered_map{}; + auto up = std::vector{}; + + // Start at well level, accumulate contributions towards root of + // group tree (FIELD group). + + for (const auto& wname : schedule_.wellNames(reportStepIdx)) { + if (! (this->wellState().hasWellRates(wname) && + this->guideRate_->has(wname))) + { + continue; + } + + const auto& well = schedule_.getWell(wname, reportStepIdx); + const auto& parent = well.groupName(); + + if (parent == "FIELD") { + // Well parented directly to "FIELD". Inadvisable and + // unexpected, but nothing to do about that here. Just skip + // this guide rate contribution. + continue; + } + + auto& grval = well.isInjector() + ? gr[parent].injection + : gr[parent].production; + + grval += this->getGuideRateValues(well); + up.push_back(parent); + } + + // Propagate accumulated guide rates up towards root of group tree. + // Override accumulation if there is a GUIDERAT specification that + // applies to a group. + std::sort(up.begin(), up.end()); + auto start = 0*up.size(); + auto u = std::unique(up.begin(), up.end()); + auto nu = std::distance(up.begin(), u); + while (nu > 0) { + const auto ntot = up.size(); + + for (auto gi = 0*nu; gi < nu; ++gi) { + const auto& gname = up[start + gi]; + const auto& group = schedule_.getGroup(gname, reportStepIdx); + + if (this->guideRate_->has(gname)) { + gr[gname].production = this->getGuideRateValues(group); + } + + if (this->guideRate_->has(gname, Phase::WATER) + || this->guideRate_->has(gname, Phase::GAS)) { + gr[gname].injection = this->getGuideRateInjectionGroupValues(group); + } + + const auto parent = group.parent(); + if (parent == "FIELD") { continue; } + + gr[parent].injection += gr[gname].injection; + gr[parent].production += gr[gname].production; + up.push_back(parent); + } + + start = ntot; + + auto begin = up.begin() + ntot; + std::sort(begin, up.end()); + u = std::unique(begin, up.end()); + nu = std::distance(begin, u); + } + + return gr; +} + +void +BlackoilWellModelGeneric:: +assignGroupControl(const Group& group, + data::GroupData& gdata) const +{ + const auto& gname = group.name(); + const auto grup_type = group.getGroupType(); + auto& cgc = gdata.currentControl; + + cgc.currentProdConstraint = Group::ProductionCMode::NONE; + + cgc.currentGasInjectionConstraint = + cgc.currentWaterInjectionConstraint = Group::InjectionCMode::NONE; + + if (this->groupState().has_production_control(gname)) { + cgc.currentProdConstraint = this->groupState().production_control(gname); + } + + if ((grup_type == ::Opm::Group::GroupType::INJECTION) || + (grup_type == ::Opm::Group::GroupType::MIXED)) + { + if (this->groupState().has_injection_control(gname, Phase::WATER)) { + cgc.currentWaterInjectionConstraint = this->groupState().injection_control(gname, Phase::WATER); + } + + if (this->groupState().has_injection_control(gname, Phase::GAS)) { + cgc.currentGasInjectionConstraint = this->groupState().injection_control(gname, Phase::GAS); + } + } +} + +void +BlackoilWellModelGeneric:: +assignGroupGuideRates(const Group& group, + const std::unordered_map& groupGuideRates, + data::GroupData& gdata) const +{ + auto& prod = gdata.guideRates.production; prod.clear(); + auto& inj = gdata.guideRates.injection; inj .clear(); + + auto xgrPos = groupGuideRates.find(group.name()); + if ((xgrPos == groupGuideRates.end()) || + !this->guideRate_->has(group.name())) + { + // No guiderates defined for this group. + return; + } + + const auto& xgr = xgrPos->second; + prod = xgr.production; + inj = xgr.injection; +} + +void +BlackoilWellModelGeneric:: +assignGroupValues(const int reportStepIdx, + std::map& gvalues) const +{ + const auto groupGuideRates = + this->calculateAllGroupGuiderates(reportStepIdx); + + for (const auto& gname : schedule_.groupNames(reportStepIdx)) { + const auto& grup = schedule_.getGroup(gname, reportStepIdx); + + auto& gdata = gvalues[gname]; + this->assignGroupControl(grup, gdata); + this->assignGroupGuideRates(grup, groupGuideRates, gdata); + } +} + +void +BlackoilWellModelGeneric:: +assignNodeValues(std::map& nodevalues) const +{ + nodevalues.clear(); + for (const auto& [node, pressure] : node_pressures_) { + nodevalues.emplace(node, data::NodeData{pressure}); + } +} + +data::GroupAndNetworkValues +BlackoilWellModelGeneric:: +groupAndNetworkData(const int reportStepIdx) const +{ + auto grp_nwrk_values = data::GroupAndNetworkValues{}; + + this->assignGroupValues(reportStepIdx, grp_nwrk_values.groupData); + this->assignNodeValues(grp_nwrk_values.nodeData); + + return grp_nwrk_values; +} + +void +BlackoilWellModelGeneric:: +updateAndCommunicateGroupData(const int reportStepIdx, + const int iterationIdx) +{ + const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); + const int nupcol = schedule()[reportStepIdx].nupcol(); + + // This builds some necessary lookup structures, so it must be called + // before we copy to well_state_nupcol_. + this->wellState().updateGlobalIsGrup(comm_); + + if (iterationIdx < nupcol) { + this->updateNupcolWGState(); + } + + auto& well_state = this->wellState(); + const auto& well_state_nupcol = this->nupcolWellState(); + // the group target reduction rates needs to be update since wells may have switched to/from GRUP control + // Currently the group target reduction does not honor NUPCOL. TODO: is that true? + std::vector groupTargetReduction(numPhases(), 0.0); + WellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, phase_usage_, *guideRate_, well_state_nupcol, well_state, this->groupState(), groupTargetReduction); + std::vector groupTargetReductionInj(numPhases(), 0.0); + WellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, phase_usage_, *guideRate_, well_state_nupcol, well_state, this->groupState(), groupTargetReductionInj); + + WellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, phase_usage_, summaryState_, well_state_nupcol, well_state, this->groupState()); + WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + + WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + + // We use the rates from the previous time-step to reduce oscillations + WellGroupHelpers::updateWellRates(fieldGroup, schedule(), reportStepIdx, this->prevWellState(), well_state); + + // Set ALQ for off-process wells to zero + for (const auto& wname : schedule().wellNames(reportStepIdx)) { + const bool is_producer = schedule().getWell(wname, reportStepIdx).isProducer(); + const bool not_on_this_process = well_state.wellMap().count(wname) == 0; + if (is_producer && not_on_this_process) { + well_state.setALQ(wname, 0.0); + } + } + + well_state.communicateGroupRates(comm_); + this->groupState().communicate_rates(comm_); + // compute wsolvent fraction for REIN wells + updateWsolvent(fieldGroup, reportStepIdx, well_state_nupcol); +} + +} diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp new file mode 100644 index 000000000..656590eb2 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -0,0 +1,355 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_BLACKOILWELLMODEL_GENERIC_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_GENERIC_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace Opm { + +namespace data { +struct GroupData; +struct GroupGuideRates; +struct GroupAndNetworkValues; +struct NodeData; +} + +class DeferredLogger; +class EclipseState; +class Group; +class RestartValue; +class Schedule; +class WellState; + +/// Class for handling the blackoil well model. +class BlackoilWellModelGeneric +{ +public: + // --------- Types --------- + using Comm = Dune::CollectiveCommunication; + + BlackoilWellModelGeneric(const Schedule& schedule, + const SummaryState& summaryState, + const EclipseState& eclState, + const PhaseUsage& phase_usage, + const Comm& comm); + + virtual ~BlackoilWellModelGeneric() = default; + + int numLocalWells() const; + int numPhases() const; + + /// return true if wells are available in the reservoir + bool wellsActive() const; + bool hasWell(const std::string& wname); + /// return true if wells are available on this process + bool localWellsActive() const; + // whether there exists any multisegment well open on this process + bool anyMSWellOpenLocal() const; + + const Well& getWellEcl(const std::string& well_name) const; + std::vector getLocalWells(const int timeStepIdx) const; + const Schedule& schedule() const { return schedule_; } + const PhaseUsage& phaseUsage() const { return phase_usage_; } + const GroupState& groupState() const { return this->active_wgstate_.group_state; } + + /* + Immutable version of the currently active wellstate. + */ + const WellState& wellState() const + { + return this->active_wgstate_.well_state; + } + + /* + Mutable version of the currently active wellstate. + */ + WellState& wellState() + { + return this->active_wgstate_.well_state; + } + + double wellPI(const int well_index) const; + double wellPI(const std::string& well_name) const; + + void updateEclWells(const int timeStepIdx, + const std::unordered_set& wells); + + + void loadRestartData(const data::Wells& rst_wells, + const data::GroupAndNetworkValues& grpNwrkValues, + const PhaseUsage& phases, + const bool handle_ms_well, + WellState& well_state); + + void initFromRestartFile(const RestartValue& restartValues, + const size_t numCells, + bool handle_ms_well); + + void setWellsActive(const bool wells_active); + + /* + Will assign the internal member last_valid_well_state_ to the + current value of the this->active_well_state_. The state stored + with storeWellState() can then subsequently be recovered with the + resetWellState() method. + */ + void commitWGState() + { + this->last_valid_wgstate_ = this->active_wgstate_; + } + + data::GroupAndNetworkValues groupAndNetworkData(const int reportStepIdx) const; + +protected: + GroupState& groupState() { return this->active_wgstate_.group_state; } + + /* + The dynamic state of the well model is maintained with an instance + of the WellState class. Currently we have + three different wellstate instances: + + 1. The currently active wellstate is in the active_well_state_ + member. That is the state which is mutated by the simulator. + + 2. In the case timestep fails to converge and we must go back and + try again with a smaller timestep we need to recover the last + valid wellstate. This is maintained with the + last_valid_well_state_ member and the functions + commitWellState() and resetWellState(). + + 3. For the NUPCOL functionality we should either use the + currently active wellstate or a wellstate frozen at max + nupcol iterations. This is handled with the member + nupcol_well_state_ and the initNupcolWellState() function. + */ + + + /* + Will return the last good wellstate. This is typcially used when + initializing a new report step where the Schedule object might + have introduced new wells. The wellstate returned by + prevWellState() must have been stored with the commitWellState() + function first. + */ + const WellState& prevWellState() const + { + return this->last_valid_wgstate_.well_state; + } + + const WGState& prevWGState() const + { + return this->last_valid_wgstate_; + } + /* + Will return the currently active nupcolWellState; must initialize + the internal nupcol wellstate with initNupcolWellState() first. + */ + const WellState& nupcolWellState() const + { + return this->nupcol_wgstate_.well_state; + } + + /* + Will store a copy of the input argument well_state in the + last_valid_well_state_ member, that state can then be recovered + with a subsequent call to resetWellState(). + */ + void commitWGState(WGState wgstate) + { + this->last_valid_wgstate_ = std::move(wgstate); + } + + /* + Will update the internal variable active_well_state_ to whatever + was stored in the last_valid_well_state_ member. This function + works in pair with commitWellState() which should be called first. + */ + void resetWGState() + { + this->active_wgstate_ = this->last_valid_wgstate_; + } + + /* + Will store the current active wellstate in the nupcol_well_state_ + member. This can then be subsequently retrieved with accessor + nupcolWellState(). + */ + void updateNupcolWGState() + { + this->nupcol_wgstate_ = this->active_wgstate_; + } + + /// \brief Create the parallel well information + /// \param localWells The local wells from ECL schedule + std::vector createLocalParallelWellInfo(const std::vector& wells); + + void initializeWellProdIndCalculators(); + void initializeWellPerfData(); + + bool wasDynamicallyShutThisTimeStep(const int well_index) const; + + void updateWsolvent(const Group& group, + const int reportStepIdx, + const WellState& wellState); + void setWsolvent(const Group& group, + const int reportStepIdx, + double wsolvent); + virtual void setWellWsolvent(const std::string& name, double wsolvent) = 0; + virtual void calcRates(const int fipnum, + const int pvtreg, + std::vector& resv_coeff) = 0; + + data::GuideRateValue getGuideRateValues(const Group& group) const; + data::GuideRateValue getGuideRateValues(const Well& well) const; + data::GuideRateValue getGuideRateInjectionGroupValues(const Group& group) const; + void getGuideRateValues(const GuideRate::RateVector& qs, + const bool is_inj, + const std::string& wgname, + data::GuideRateValue& grval) const; + + void assignWellGuideRates(data::Wells& wsrpt) const; + void assignShutConnections(data::Wells& wsrpt, + const int reportStepIndex) const; + void assignGroupControl(const Group& group, + data::GroupData& gdata) const; + void assignGroupGuideRates(const Group& group, + const std::unordered_map& groupGuideRates, + data::GroupData& gdata) const; + void assignGroupValues(const int reportStepIdx, + std::map& gvalues) const; + void assignNodeValues(std::map& nodevalues) const; + + std::unordered_map + calculateAllGroupGuiderates(const int reportStepIdx) const; + + bool checkGroupConstraints(const Group& group, + const int reportStepIdx, + DeferredLogger& deferred_logger) const; + + Group::InjectionCMode checkGroupInjectionConstraints(const Group& group, + const int reportStepIdx, + const Phase& phase) const; + Group::ProductionCMode checkGroupProductionConstraints(const Group& group, + const int reportStepIdx, + DeferredLogger& deferred_logger) const; + + void checkGconsaleLimits(const Group& group, + WellState& well_state, + const int reportStepIdx, + DeferredLogger& deferred_logger); + + void checkGroupHigherConstraints(const Group& group, + DeferredLogger& deferred_logger, + const int reportStepIdx, + std::set& switched_groups); + + void updateGroupIndividualControl(const Group& group, + DeferredLogger& deferred_logger, + const int reportStepIdx, + std::set& switched_groups); + + void updateGroupIndividualControls(DeferredLogger& deferred_logger, + std::set& switched_groups, + const int reportStepIdx, + const int iterationIdx); + + void updateGroupHigherControls(DeferredLogger& deferred_logger, + const int reportStepIdx, + std::set& switched_groups); + + void actionOnBrokenConstraints(const Group& group, + const Group::ExceedAction& exceed_action, + const Group::ProductionCMode& newControl, + DeferredLogger& deferred_logger); + void actionOnBrokenConstraints(const Group& group, + const Group::InjectionCMode& newControl, + const Phase& controlPhase, + DeferredLogger& deferred_logger); + + void updateAndCommunicateGroupData(const int reportStepIdx, + const int iterationIdx); + + const Schedule& schedule_; + const SummaryState& summaryState_; + const EclipseState& eclState_; + const Comm& comm_; + + PhaseUsage phase_usage_; + bool terminal_output_{false}; + bool wells_active_{false}; + bool initial_step_{}; + + std::vector wells_ecl_; + std::vector> well_perf_data_; + std::function not_on_process_{}; + + std::vector parallel_well_info_; + std::vector local_parallel_well_info_; + + std::vector prod_index_calc_; + + // Map from logically cartesian cell indices to compressed ones. + // Cells not in the interior are not mapped. This deactivates + // these for distributed wells and makes the distribution non-overlapping. + std::vector cartesian_to_compressed_; + + std::vector pvt_region_idx_; + + mutable std::unordered_set closed_this_step_; + + WellTestState wellTestState_{}; + std::unique_ptr guideRate_; + std::map node_pressures_; // Storing network pressures for output. + + /* + The various wellState members should be accessed and modified + through the accessor functions wellState(), prevWellState(), + commitWellState(), resetWellState(), nupcolWellState() and + updateNupcolWellState(). + */ + WGState active_wgstate_; + WGState last_valid_wgstate_; + WGState nupcol_wgstate_; +}; + + +} // namespace Opm + +#endif diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index d87d43037..8864d18fd 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -34,17 +34,15 @@ namespace Opm { template BlackoilWellModel:: BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage) - : ebosSimulator_(ebosSimulator) - , terminal_output_((ebosSimulator.gridView().comm().rank() == 0) && - EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput)) - , phase_usage_(phase_usage) - , active_wgstate_(phase_usage) - , last_valid_wgstate_(phase_usage) - , nupcol_wgstate_(phase_usage) + : BlackoilWellModelGeneric(ebosSimulator.vanguard().schedule(), + ebosSimulator.vanguard().summaryState(), + ebosSimulator.vanguard().eclState(), + phase_usage, + ebosSimulator.gridView().comm()) + , ebosSimulator_(ebosSimulator) { - // Create the guide rate container. - this->guideRate_ = - std::make_unique(ebosSimulator_.vanguard().schedule()); + terminal_output_ = ((ebosSimulator.gridView().comm().rank() == 0) && + EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput)); local_num_cells_ = ebosSimulator_.gridView().size(0); @@ -63,21 +61,6 @@ namespace Opm { parallel_wells.end()); } - const auto numProcs = ebosSimulator.gridView().comm().size(); - this->not_on_process_ = [this, numProcs](const Well& well) { - if (numProcs == decltype(numProcs){1}) - return false; - - // Recall: false indicates NOT active! - const auto value = std::make_pair(well.name(), true); - auto candidate = std::lower_bound(this->parallel_well_info_.begin(), - this->parallel_well_info_.end(), - value); - - return (candidate == this->parallel_well_info_.end()) - || (*candidate != value); - }; - this->alternative_well_rate_init_ = EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit); } @@ -226,34 +209,6 @@ namespace Opm { } - template - std::vector< Well > - BlackoilWellModel:: - getLocalWells(const int timeStepIdx) const - { - auto w = schedule().getWells(timeStepIdx); - w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end()); - return w; - } - - template - std::vector< ParallelWellInfo* > - BlackoilWellModel::createLocalParallelWellInfo(const std::vector& wells) - { - std::vector< ParallelWellInfo* > local_parallel_well_info; - local_parallel_well_info.reserve(wells.size()); - for (const auto& well : wells) - { - auto wellPair = std::make_pair(well.name(), true); - auto pwell = std::lower_bound(parallel_well_info_.begin(), - parallel_well_info_.end(), - wellPair); - assert(pwell != parallel_well_info_.end() && - *pwell == wellPair); - local_parallel_well_info.push_back(&(*pwell)); - } - return local_parallel_well_info; - } template void @@ -319,9 +274,10 @@ namespace Opm { DeferredLogger local_deferredLogger; this->resetWGState(); - updateAndCommunicateGroupData(); - this->wellState().gliftTimeStepInit(); const int reportStepIdx = ebosSimulator_.episodeIndex(); + updateAndCommunicateGroupData(reportStepIdx, + ebosSimulator_.model().newtonMethod().numIterations()); + this->wellState().gliftTimeStepInit(); const double simulationTime = ebosSimulator_.time(); std::string exc_msg; auto exc_type = ExceptionType::NONE; @@ -566,7 +522,8 @@ namespace Opm { // check group sales limits at the end of the timestep const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); - checkGconsaleLimits(fieldGroup, this->wellState(), local_deferredLogger); + checkGconsaleLimits(fieldGroup, this->wellState(), + ebosSimulator_.episodeIndex(), local_deferredLogger); this->calculateProductivityIndexValues(local_deferredLogger); @@ -617,119 +574,6 @@ namespace Opm { return nullptr; } - template - void - BlackoilWellModel:: - initFromRestartFile(const RestartValue& restartValues) - { - // The restart step value is used to identify wells present at the given - // time step. Wells that are added at the same time step as RESTART is initiated - // will not be present in a restart file. Use the previous time step to retrieve - // wells that have information written to the restart file. - const int report_step = std::max(eclState().getInitConfig().getRestartStep() - 1, 0); - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - // wells_ecl_ should only contain wells on this processor. - wells_ecl_ = getLocalWells(report_step); - local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_); - - this->initializeWellProdIndCalculators(); - initializeWellPerfData(); - - const int nw = wells_ecl_.size(); - if (nw > 0) { - const auto phaseUsage = phaseUsageFromDeck(eclState()); - const size_t numCells = UgGridHelpers::numCells(grid()); - const bool handle_ms_well = (param_.use_multisegment_well_ && anyMSWellOpenLocal()); - this->wellState().resize(wells_ecl_, local_parallel_well_info_, schedule(), handle_ms_well, numCells, well_perf_data_, summaryState); // Resize for restart step - loadRestartData(restartValues.wells, restartValues.grp_nwrk, phaseUsage, handle_ms_well, this->wellState()); - } - - this->commitWGState(); - initial_step_ = false; - } - - - - - - template - void - BlackoilWellModel:: - initializeWellProdIndCalculators() - { - this->prod_index_calc_.clear(); - this->prod_index_calc_.reserve(this->wells_ecl_.size()); - for (const auto& well : this->wells_ecl_) { - this->prod_index_calc_.emplace_back(well); - } - } - - - - - - template - void - BlackoilWellModel:: - initializeWellPerfData() - { - well_perf_data_.resize(wells_ecl_.size()); - int well_index = 0; - for (const auto& well : wells_ecl_) { - int completion_index = 0; - // INVALID_ECL_INDEX marks no above perf available - int completion_index_above = ParallelWellInfo::INVALID_ECL_INDEX; - well_perf_data_[well_index].clear(); - well_perf_data_[well_index].reserve(well.getConnections().size()); - CheckDistributedWellConnections checker(well, *local_parallel_well_info_[well_index]); - bool hasFirstPerforation = false; - bool firstOpenCompletion = true; - auto& parallelWellInfo = *local_parallel_well_info_[well_index]; - parallelWellInfo.beginReset(); - - for (const auto& completion : well.getConnections()) { - const int active_index = - cartesian_to_compressed_[completion.global_index()]; - if (completion.state() == Connection::State::OPEN) { - if (active_index >= 0) { - if (firstOpenCompletion) - { - hasFirstPerforation = true; - } - checker.connectionFound(completion_index); - PerforationData pd; - pd.cell_index = active_index; - pd.connection_transmissibility_factor = completion.CF(); - pd.satnum_id = completion.satTableId(); - pd.ecl_index = completion_index; - well_perf_data_[well_index].push_back(pd); - parallelWellInfo.pushBackEclIndex(completion_index_above, - completion_index); - } - firstOpenCompletion = false; - // Next time this index is the one above as each open completion is - // is stored somehwere. - completion_index_above = completion_index; - } else { - checker.connectionFound(completion_index); - if (completion.state() != Connection::State::SHUT) { - OPM_THROW(std::runtime_error, - "Completion state: " << Connection::State2String(completion.state()) << " not handled"); - } - } - // Note: we rely on the connections being filtered! I.e. there are only connections - // to active cells in the global grid. - ++completion_index; - } - parallelWellInfo.endReset(); - checker.checkAllConnectionsFound(); - parallelWellInfo.communicateFirstPerforation(hasFirstPerforation); - ++well_index; - } - } - - - template @@ -1271,42 +1115,6 @@ namespace Opm { - template - bool - BlackoilWellModel:: - wellsActive() const - { - return wells_active_; - } - - - - - - template - void - BlackoilWellModel:: - setWellsActive(const bool wells_active) - { - wells_active_ = wells_active; - } - - - - - - template - bool - BlackoilWellModel:: - localWellsActive() const - { - return numLocalWells() > 0; - } - - - - - template void BlackoilWellModel:: @@ -1357,7 +1165,9 @@ namespace Opm { if (checkGroupConvergence) { const int reportStepIdx = ebosSimulator_.episodeIndex(); const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); - bool violated = checkGroupConstraints(fieldGroup, global_deferredLogger); + bool violated = checkGroupConstraints(fieldGroup, + ebosSimulator_.episodeIndex(), + global_deferredLogger); report.setGroupConverged(!violated); } return report; @@ -1392,7 +1202,9 @@ namespace Opm { // For no well active globally we simply return. if( !wellsActive() ) return ; - updateAndCommunicateGroupData(); + const int episodeIdx = ebosSimulator_.episodeIndex(); + const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations(); + updateAndCommunicateGroupData(episodeIdx, iterationIdx); updateNetworkPressures(); @@ -1401,12 +1213,15 @@ namespace Opm { if (checkGroupControls) { // Check group individual constraints. - updateGroupIndividualControls(deferred_logger, switched_groups); + updateGroupIndividualControls(deferred_logger, switched_groups, + episodeIdx, iterationIdx); // Check group's constraints from higher levels. - updateGroupHigherControls(deferred_logger, switched_groups); + updateGroupHigherControls(deferred_logger, + episodeIdx, + switched_groups); - updateAndCommunicateGroupData(); + updateAndCommunicateGroupData(episodeIdx, iterationIdx); // Check wells' group constraints and communicate. for (const auto& well : well_container_) { @@ -1416,7 +1231,7 @@ namespace Opm { switched_wells.insert(well->name()); } } - updateAndCommunicateGroupData(); + updateAndCommunicateGroupData(episodeIdx, iterationIdx); } // Check individual well constraints and communicate. @@ -1427,8 +1242,7 @@ namespace Opm { const auto mode = WellInterface::IndividualOrGroup::Individual; well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger); } - updateAndCommunicateGroupData(); - + updateAndCommunicateGroupData(episodeIdx, iterationIdx); } @@ -1466,64 +1280,6 @@ namespace Opm { - template - void - BlackoilWellModel:: - updateAndCommunicateGroupData() - { - const int reportStepIdx = ebosSimulator_.episodeIndex(); - const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); - const int nupcol = schedule()[reportStepIdx].nupcol(); - const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations(); - - // This builds some necessary lookup structures, so it must be called - // before we copy to well_state_nupcol_. - const auto& comm = ebosSimulator_.vanguard().grid().comm(); - this->wellState().updateGlobalIsGrup(comm); - - if (iterationIdx < nupcol) { - this->updateNupcolWGState(); - } - - auto& well_state = this->wellState(); - const auto& well_state_nupcol = this->nupcolWellState(); - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - // the group target reduction rates needs to be update since wells may have swicthed to/from GRUP control - // Currently the group target reduction does not honor NUPCOL. TODO: is that true? - std::vector groupTargetReduction(numPhases(), 0.0); - WellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ false, phase_usage_, *guideRate_, well_state_nupcol, well_state, this->groupState(), groupTargetReduction); - std::vector groupTargetReductionInj(numPhases(), 0.0); - WellGroupHelpers::updateGroupTargetReduction(fieldGroup, schedule(), reportStepIdx, /*isInjector*/ true, phase_usage_, *guideRate_, well_state_nupcol, well_state, this->groupState(), groupTargetReductionInj); - - WellGroupHelpers::updateREINForGroups(fieldGroup, schedule(), reportStepIdx, phase_usage_, summaryState, well_state_nupcol, well_state, this->groupState()); - WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); - - WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); - WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); - - // We use the rates from the privious time-step to reduce oscilations - WellGroupHelpers::updateWellRates(fieldGroup, schedule(), reportStepIdx, this->prevWellState(), well_state); - - // Set ALQ for off-process wells to zero - for (const auto& wname : schedule().wellNames(reportStepIdx)) { - const bool is_producer = schedule().getWell(wname, reportStepIdx).isProducer(); - const bool not_on_this_process = well_state.wellMap().count(wname) == 0; - if (is_producer && not_on_this_process) { - well_state.setALQ(wname, 0.0); - } - } - - well_state.communicateGroupRates(comm); - this->groupState().communicate_rates(comm); - // compute wsolvent fraction for REIN wells - updateWsolvent(fieldGroup, schedule(), reportStepIdx, well_state_nupcol); - - } - - - - - template void BlackoilWellModel:: @@ -1888,20 +1644,6 @@ namespace Opm { return numComp; } - template - int - BlackoilWellModel:: numLocalWells() const - { - return wells_ecl_.size(); - } - - template - int - BlackoilWellModel::numPhases() const - { - return phase_usage_.num_phases; - } - template void BlackoilWellModel::extractLegacyDepth_() @@ -1938,164 +1680,6 @@ namespace Opm { } - template - bool - BlackoilWellModel:: - hasWell(const std::string& wname) { - auto iter = std::find_if( this->wells_ecl_.begin(), this->wells_ecl_.end(), [&wname](const Well& well) { return well.name() == wname; }); - return (iter != this->wells_ecl_.end()); - } - - - // convert well data from opm-common to well state from opm-core - template - void - BlackoilWellModel:: - loadRestartData( const data::Wells& rst_wells, - const data::GroupAndNetworkValues& grpNwrkValues, - const PhaseUsage& phases, - const bool handle_ms_well, - WellState& well_state) - { - using GPMode = Group::ProductionCMode; - using GIMode = Group::InjectionCMode; - - using rt = data::Rates::opt; - const auto np = phases.num_phases; - - std::vector< rt > phs( np ); - if( phases.phase_used[BlackoilPhases::Aqua] ) { - phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat; - } - - if( phases.phase_used[BlackoilPhases::Liquid] ) { - phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil; - } - - if( phases.phase_used[BlackoilPhases::Vapour] ) { - phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas; - } - - for( const auto& wm : well_state.wellMap() ) { - const auto well_index = wm.second[ 0 ]; - const auto& rst_well = rst_wells.at( wm.first ); - well_state.update_bhp( well_index, rst_well.bhp); - well_state.update_temperature( well_index, rst_well.temperature); - - if (rst_well.current_control.isProducer) { - well_state.currentProductionControl( well_index, rst_well.current_control.prod); - } - else { - well_state.currentInjectionControl( well_index, rst_well.current_control.inj); - } - - for( size_t i = 0; i < phs.size(); ++i ) { - assert( rst_well.rates.has( phs[ i ] ) ); - well_state.wellRates(well_index)[ i ] = rst_well.rates.get( phs[ i ] ); - } - - auto& perf_pressure = well_state.perfPress(well_index); - auto& perf_rates = well_state.perfRates(well_index); - auto& perf_phase_rates = well_state.perfPhaseRates(well_index); - const auto& perf_data = this->well_perf_data_[well_index]; - - for (std::size_t perf_index = 0; perf_index < perf_data.size(); perf_index++) { - const auto& pd = perf_data[perf_index]; - const auto& rst_connection = rst_well.connections[pd.ecl_index]; - perf_pressure[perf_index] = rst_connection.pressure; - perf_rates[perf_index] = rst_connection.reservoir_rate; - for (int phase_index = 0; phase_index < np; ++phase_index) - perf_phase_rates[perf_index*np + phase_index] = rst_connection.rates.get(phs[phase_index]); - } - - if (handle_ms_well && !rst_well.segments.empty()) { - // we need the well_ecl_ information - const std::string& well_name = wm.first; - const Well& well_ecl = getWellEcl(well_name); - - const WellSegments& segment_set = well_ecl.getSegments(); - - const auto& rst_segments = rst_well.segments; - - // \Note: eventually we need to hanlde the situations that some segments are shut - assert(0u + segment_set.size() == rst_segments.size()); - - auto& segments = well_state.segments(well_index); - auto& segment_pressure = segments.pressure; - auto& segment_rates = segments.rates; - for (const auto& rst_segment : rst_segments) { - const int segment_index = segment_set.segmentNumberToIndex(rst_segment.first); - - // recovering segment rates and pressure from the restart values - const auto pres_idx = data::SegmentPressures::Value::Pressure; - segment_pressure[segment_index] = rst_segment.second.pressures[pres_idx]; - - const auto& rst_segment_rates = rst_segment.second.rates; - for (int p = 0; p < np; ++p) { - segment_rates[segment_index * np + p] = rst_segment_rates.get(phs[p]); - } - } - } - } - - for (const auto& [group, value] : grpNwrkValues.groupData) { - const auto cpc = value.currentControl.currentProdConstraint; - const auto cgi = value.currentControl.currentGasInjectionConstraint; - const auto cwi = value.currentControl.currentWaterInjectionConstraint; - - if (cpc != GPMode::NONE) { - this->groupState().production_control(group, cpc); - } - - if (cgi != GIMode::NONE) { - this->groupState().injection_control(group, Phase::GAS, cgi); - } - - if (cwi != GIMode::NONE) { - this->groupState().injection_control(group, Phase::WATER, cwi); - } - } - } - - - - - - template - bool - BlackoilWellModel:: - anyMSWellOpenLocal() const - { - for (const auto& well : wells_ecl_) { - if (well.isMultiSegment()) { - return true; - } - } - return false; - } - - - - - - template - const Well& - BlackoilWellModel:: - getWellEcl(const std::string& well_name) const - { - // finding the iterator of the well in wells_ecl - auto well_ecl = std::find_if(wells_ecl_.begin(), - wells_ecl_.end(), - [&well_name](const Well& elem)->bool { - return elem.name() == well_name; - }); - - assert(well_ecl != wells_ecl_.end()); - - return *well_ecl; - } - - template typename BlackoilWellModel::WellInterfacePtr BlackoilWellModel:: @@ -2116,731 +1700,6 @@ namespace Opm { - template - void - BlackoilWellModel:: - updateGroupIndividualControls(DeferredLogger& deferred_logger, std::set& switched_groups) - { - const int reportStepIdx = ebosSimulator_.episodeIndex(); - - const int nupcol = schedule()[reportStepIdx].nupcol(); - const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations(); - // don't switch group control when iterationIdx > nupcol - // to avoid oscilations between group controls - if (iterationIdx > nupcol) - return; - - const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); - updateGroupIndividualControl(fieldGroup, deferred_logger, switched_groups); - } - - - - template - void - BlackoilWellModel:: - updateGroupIndividualControl(const Group& group, DeferredLogger& deferred_logger, std::set& switched_groups) { - - const int reportStepIdx = ebosSimulator_.episodeIndex(); - const bool skip = switched_groups.count(group.name()); - if (!skip && group.isInjectionGroup()) - { - const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; - for (Phase phase : all) { - if (!group.hasInjectionControl(phase)) { - continue; - } - Group::InjectionCMode newControl = checkGroupInjectionConstraints(group, phase); - if (newControl != Group::InjectionCMode::NONE) - { - switched_groups.insert(group.name()); - actionOnBrokenConstraints(group, newControl, phase, deferred_logger); - } - } - } - if (!skip && group.isProductionGroup()) { - Group::ProductionCMode newControl = checkGroupProductionConstraints(group, deferred_logger); - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - const auto controls = group.productionControls(summaryState); - if (newControl != Group::ProductionCMode::NONE) - { - switched_groups.insert(group.name()); - actionOnBrokenConstraints(group, controls.exceed_action, newControl, deferred_logger); - } - } - - // call recursively down the group hiearchy - for (const std::string& groupName : group.groups()) { - updateGroupIndividualControl( schedule().getGroup(groupName, reportStepIdx), deferred_logger, switched_groups); - } - } - - template - bool - BlackoilWellModel:: - checkGroupConstraints(const Group& group, DeferredLogger& deferred_logger) const { - - const int reportStepIdx = ebosSimulator_.episodeIndex(); - if (group.isInjectionGroup()) { - const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; - for (Phase phase : all) { - if (!group.hasInjectionControl(phase)) { - continue; - } - Group::InjectionCMode newControl = checkGroupInjectionConstraints(group, phase); - if (newControl != Group::InjectionCMode::NONE) { - return true; - } - } - } - if (group.isProductionGroup()) { - Group::ProductionCMode newControl = checkGroupProductionConstraints(group, deferred_logger); - if (newControl != Group::ProductionCMode::NONE) - { - return true; - } - } - - // call recursively down the group hiearchy - bool violated = false; - for (const std::string& groupName : group.groups()) { - violated = violated || checkGroupConstraints( schedule().getGroup(groupName, reportStepIdx), deferred_logger); - } - return violated; - } - - - - template - Group::ProductionCMode - BlackoilWellModel:: - checkGroupProductionConstraints(const Group& group, DeferredLogger& deferred_logger) const { - - const int reportStepIdx = ebosSimulator_.episodeIndex(); - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - const auto& comm = ebosSimulator_.vanguard().grid().comm(); - const auto& well_state = this->wellState(); - - const auto controls = group.productionControls(summaryState); - const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name()); - - if (group.has_control(Group::ProductionCMode::ORAT)) - { - if (currentControl != Group::ProductionCMode::ORAT) - { - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false); - - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.oil_target < current_rate ) { - return Group::ProductionCMode::ORAT; - } - } - } - - if (group.has_control(Group::ProductionCMode::WRAT)) - { - if (currentControl != Group::ProductionCMode::WRAT) - { - - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false); - - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.water_target < current_rate ) { - return Group::ProductionCMode::WRAT; - } - } - } - if (group.has_control(Group::ProductionCMode::GRAT)) - { - if (currentControl != Group::ProductionCMode::GRAT) - { - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false); - - // sum over all nodes - current_rate = comm.sum(current_rate); - if (controls.gas_target < current_rate ) { - return Group::ProductionCMode::GRAT; - } - } - } - if (group.has_control(Group::ProductionCMode::LRAT)) - { - if (currentControl != Group::ProductionCMode::LRAT) - { - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false); - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false); - - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.liquid_target < current_rate ) { - return Group::ProductionCMode::LRAT; - } - } - } - - if (group.has_control(Group::ProductionCMode::CRAT)) - { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "CRAT control for production groups not implemented" , deferred_logger); - } - if (group.has_control(Group::ProductionCMode::RESV)) - { - if (currentControl != Group::ProductionCMode::RESV) - { - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true); - current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true); - current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true); - - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.resv_target < current_rate ) { - return Group::ProductionCMode::RESV; - } - } - - } - if (group.has_control(Group::ProductionCMode::PRBL)) - { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "PRBL control for production groups not implemented", deferred_logger); - } - return Group::ProductionCMode::NONE; - } - - - template - Group::InjectionCMode - BlackoilWellModel:: - checkGroupInjectionConstraints(const Group& group, const Phase& phase) const { - - const int reportStepIdx = ebosSimulator_.episodeIndex(); - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - const auto& comm = ebosSimulator_.vanguard().grid().comm(); - const auto& well_state = this->wellState(); - - int phasePos; - if (phase == Phase::GAS && phase_usage_.phase_used[BlackoilPhases::Vapour] ) - phasePos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; - else if (phase == Phase::OIL && phase_usage_.phase_used[BlackoilPhases::Liquid]) - phasePos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; - else if (phase == Phase::WATER && phase_usage_.phase_used[BlackoilPhases::Aqua] ) - phasePos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - else - OPM_THROW(std::runtime_error, "Unknown phase" ); - - const auto& controls = group.injectionControls(phase, summaryState); - auto currentControl = this->groupState().injection_control(group.name(), phase); - - if (controls.has_control(Group::InjectionCMode::RATE)) - { - if (currentControl != Group::InjectionCMode::RATE) - { - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true); - - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.surface_max_rate < current_rate) { - return Group::InjectionCMode::RATE; - } - } - } - if (controls.has_control(Group::InjectionCMode::RESV)) - { - if (currentControl != Group::InjectionCMode::RESV) - { - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true); - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.resv_max_rate < current_rate) { - return Group::InjectionCMode::RESV; - } - } - } - if (controls.has_control(Group::InjectionCMode::REIN)) - { - if (currentControl != Group::InjectionCMode::REIN) - { - double production_Rate = 0.0; - const Group& groupRein = schedule().getGroup(controls.reinj_group, reportStepIdx); - production_Rate += WellGroupHelpers::sumWellRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false); - - // sum over all nodes - production_Rate = comm.sum(production_Rate); - - double current_rate = 0.0; - current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true); - - // sum over all nodes - current_rate = comm.sum(current_rate); - - if (controls.target_reinj_fraction*production_Rate < current_rate) { - return Group::InjectionCMode::REIN; - } - } - } - if (controls.has_control(Group::InjectionCMode::VREP)) - { - if (currentControl != Group::InjectionCMode::VREP) - { - double voidage_rate = 0.0; - const Group& groupVoidage = schedule().getGroup(controls.voidage_group, reportStepIdx); - voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false); - voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false); - voidage_rate += WellGroupHelpers::sumWellResRates(groupVoidage, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false); - - // sum over all nodes - voidage_rate = comm.sum(voidage_rate); - - double total_rate = 0.0; - total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], true); - total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], true); - total_rate += WellGroupHelpers::sumWellResRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], true); - - // sum over all nodes - total_rate = comm.sum(total_rate); - - if (controls.target_void_fraction*voidage_rate < total_rate) { - return Group::InjectionCMode::VREP; - } - } - } - return Group::InjectionCMode::NONE; - } - - template - void - BlackoilWellModel:: - checkGconsaleLimits(const Group& group, WellState& well_state, DeferredLogger& deferred_logger) - { - const int reportStepIdx = ebosSimulator_.episodeIndex(); - // call recursively down the group hiearchy - for (const std::string& groupName : group.groups()) { - checkGconsaleLimits( schedule().getGroup(groupName, reportStepIdx), well_state, deferred_logger); - } - - // only for groups with gas injection controls - if (!group.hasInjectionControl(Phase::GAS)) { - return; - } - - // check if gconsale is used for this group - if (!schedule()[reportStepIdx].gconsale().has(group.name())) - return; - - std::ostringstream ss; - - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - const auto& comm = ebosSimulator_.vanguard().grid().comm(); - - const auto& gconsale = schedule()[reportStepIdx].gconsale().get(group.name(), summaryState); - const Group::ProductionCMode& oldProductionControl = this->groupState().production_control(group.name()); - - - int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; - double production_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/false); - double injection_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/true); - - // sum over all nodes - injection_rate = comm.sum(injection_rate); - production_rate = comm.sum(production_rate); - - double sales_rate = production_rate - injection_rate; - double production_target = gconsale.sales_target + injection_rate; - - // add import rate and substract consumption rate for group for gas - if (schedule()[reportStepIdx].gconsump().has(group.name())) { - const auto& gconsump = schedule()[reportStepIdx].gconsump().get(group.name(), summaryState); - if (phase_usage_.phase_used[BlackoilPhases::Vapour]) { - sales_rate += gconsump.import_rate; - sales_rate -= gconsump.consumption_rate; - production_target -= gconsump.import_rate; - production_target += gconsump.consumption_rate; - } - } - if (sales_rate > gconsale.max_sales_rate) { - - switch(gconsale.max_proc) { - case GConSale::MaxProcedure::NONE: { - if (oldProductionControl != Group::ProductionCMode::GRAT && oldProductionControl != Group::ProductionCMode::NONE) { - ss << "Group sales exceed maximum limit, but the action is NONE for " + group.name() + ". Nothing happens"; - } - break; - } - case GConSale::MaxProcedure::CON: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit CON not implemented", deferred_logger); - break; - } - case GConSale::MaxProcedure::CON_P: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit CON_P not implemented", deferred_logger); - break; - } - case GConSale::MaxProcedure::WELL: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit WELL not implemented", deferred_logger); - break; - } - case GConSale::MaxProcedure::PLUG: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit PLUG not implemented", deferred_logger); - break; - } - case GConSale::MaxProcedure::MAXR: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit MAXR not implemented", deferred_logger); - break; - } - case GConSale::MaxProcedure::END: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GCONSALE exceed limit END not implemented", deferred_logger); - break; - } - case GConSale::MaxProcedure::RATE: { - this->groupState().production_control(group.name(), Group::ProductionCMode::GRAT); - ss << "Maximum GCONSALE limit violated for " << group.name() << ". The group is switched from "; - ss << Group::ProductionCMode2String(oldProductionControl) << " to " << Group::ProductionCMode2String(Group::ProductionCMode::GRAT); - ss << " and limited by the maximum sales rate after consumption and import are considered" ; - this->groupState().update_grat_sales_target(group.name(), production_target); - break; - } - default: - throw("Invalid procedure for maximum rate limit selected for group" + group.name()); - } - } - if (sales_rate < gconsale.min_sales_rate) { - const Group::ProductionCMode& currentProductionControl = this->groupState().production_control(group.name()); - if ( currentProductionControl == Group::ProductionCMode::GRAT ) { - ss << "Group " + group.name() + " has sale rate less then minimum permitted value and is under GRAT control. \n"; - ss << "The GRAT is increased to meet the sales minimum rate. \n"; - this->groupState().update_grat_sales_target(group.name(), production_target); - //} else if () {//TODO add action for WGASPROD - //} else if () {//TODO add action for drilling queue - } else { - ss << "Group " + group.name() + " has sale rate less then minimum permitted value but cannot increase the group production rate \n"; - ss << "or adjust gas production using WGASPROD or drill new wells to meet the sales target. \n"; - ss << "Note that WGASPROD and drilling queues are not implemented in Flow. No action is taken. \n "; - } - } - if (gconsale.sales_target < 0.0) { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + " has sale rate target less then zero. Not implemented in Flow" , deferred_logger); - } - - auto cc = Dune::MPIHelper::getCollectiveCommunication(); - if (!ss.str().empty() && cc.rank() == 0) - deferred_logger.info(ss.str()); - - } - - - template - void - BlackoilWellModel:: - actionOnBrokenConstraints(const Group& group, const Group::ExceedAction& exceed_action, const Group::ProductionCMode& newControl, DeferredLogger& deferred_logger) { - - const Group::ProductionCMode oldControl = this->groupState().production_control(group.name()); - - std::ostringstream ss; - - switch(exceed_action) { - case Group::ExceedAction::NONE: { - if (oldControl != newControl && oldControl != Group::ProductionCMode::NONE) { - ss << "Group production exceed action is NONE for group " + group.name() + ". Nothing happens."; - } - break; - } - case Group::ExceedAction::CON: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit CON not implemented", deferred_logger); - break; - } - case Group::ExceedAction::CON_PLUS: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit CON_PLUS not implemented", deferred_logger); - break; - } - case Group::ExceedAction::WELL: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit WELL not implemented", deferred_logger); - break; - } - case Group::ExceedAction::PLUG: { - OPM_DEFLOG_THROW(std::runtime_error, "Group " + group.name() + "GroupProductionExceedLimit PLUG not implemented", deferred_logger); - break; - } - case Group::ExceedAction::RATE: { - if (oldControl != newControl) { - this->groupState().production_control(group.name(), newControl); - ss << "Switching production control mode for group "<< group.name() - << " from " << Group::ProductionCMode2String(oldControl) - << " to " << Group::ProductionCMode2String(newControl); - } - break; - } - default: - throw("Invalid procedure for maximum rate limit selected for group" + group.name()); - } - - auto cc = Dune::MPIHelper::getCollectiveCommunication(); - if (!ss.str().empty() && cc.rank() == 0) - deferred_logger.info(ss.str()); - - - } - - - template - void - BlackoilWellModel:: - actionOnBrokenConstraints(const Group& group, const Group::InjectionCMode& newControl, const Phase& controlPhase, DeferredLogger& deferred_logger) { - auto oldControl = this->groupState().injection_control(group.name(), controlPhase); - - std::ostringstream ss; - if (oldControl != newControl) { - const std::string from = Group::InjectionCMode2String(oldControl); - ss << "Switching injection control mode for group "<< group.name() - << " from " << Group::InjectionCMode2String(oldControl) - << " to " << Group::InjectionCMode2String(newControl); - this->groupState().injection_control(group.name(), controlPhase, newControl); - } - auto cc = Dune::MPIHelper::getCollectiveCommunication(); - if (!ss.str().empty() && cc.rank() == 0) - deferred_logger.info(ss.str()); - - } - - - - - template - void - BlackoilWellModel:: - updateGroupHigherControls(DeferredLogger& deferred_logger, std::set& switched_groups) - { - const int reportStepIdx = ebosSimulator_.episodeIndex(); - const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); - checkGroupHigherConstraints(fieldGroup, deferred_logger, switched_groups); - } - - - template - void - BlackoilWellModel:: - checkGroupHigherConstraints(const Group& group, DeferredLogger& deferred_logger, std::set& switched_groups) - { - // Set up coefficients for RESV <-> surface rate conversion. - // Use the pvtRegionIdx from the top cell of the first well. - // TODO fix this! - // This is only used for converting RESV rates. - // What is the proper approach? - const auto& comm = ebosSimulator_.vanguard().grid().comm(); - const int fipnum = 0; - int pvtreg = well_perf_data_.empty() || well_perf_data_[0].empty() - ? pvt_region_idx_[0] - : pvt_region_idx_[well_perf_data_[0][0].cell_index]; - - if ( comm.size() > 1) - { - // Just like in the sequential case the pvtregion is determined - // by the first cell of the first well. What is the first well - // is decided by the order in the Schedule using Well::seqIndex() - int firstWellIndex = well_perf_data_.empty() ? - std::numeric_limits::max() : wells_ecl_[0].seqIndex(); - auto regIndexPair = std::make_pair(pvtreg, firstWellIndex); - std::vector pairs(comm.size()); - comm.allgather(®IndexPair, 1, pairs.data()); - pvtreg = std::min_element(pairs.begin(), pairs.end(), - [](const auto& p1, const auto& p2){ return p1.second < p2.second;}) - ->first; - } - - std::vector resv_coeff(phase_usage_.num_phases, 0.0); - rateConverter_->calcCoeff(fipnum, pvtreg, resv_coeff); - - const int reportStepIdx = ebosSimulator_.episodeIndex(); - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - - std::vector rates(phase_usage_.num_phases, 0.0); - - const bool skip = switched_groups.count(group.name()) || group.name() == "FIELD"; - - if (!skip && group.isInjectionGroup()) { - // Obtain rates for group. - for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) { - const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ true); - // Sum over all processes - rates[phasePos] = comm.sum(local_current_rate); - } - const Phase all[] = { Phase::WATER, Phase::OIL, Phase::GAS }; - for (Phase phase : all) { - // Check higher up only if under individual (not FLD) control. - auto currentControl = this->groupState().injection_control(group.name(), phase); - if (currentControl != Group::InjectionCMode::FLD && group.injectionGroupControlAvailable(phase)) { - const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx); - const std::pair changed = WellGroupHelpers::checkGroupConstraintsInj( - group.name(), - group.parent(), - parentGroup, - this->wellState(), - this->groupState(), - reportStepIdx, - guideRate_.get(), - rates.data(), - phase, - phase_usage_, - group.getGroupEfficiencyFactor(), - schedule(), - summaryState, - resv_coeff, - deferred_logger); - if (changed.first) { - switched_groups.insert(group.name()); - actionOnBrokenConstraints(group, Group::InjectionCMode::FLD, phase, deferred_logger); - } - } - } - } - - if (!skip && group.isProductionGroup()) { - // Obtain rates for group. - for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) { - const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ false); - // Sum over all processes - rates[phasePos] = -comm.sum(local_current_rate); - } - // Check higher up only if under individual (not FLD) control. - const Group::ProductionCMode& currentControl = this->groupState().production_control(group.name()); - if (currentControl != Group::ProductionCMode::FLD && group.productionGroupControlAvailable()) { - const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx); - const std::pair changed = WellGroupHelpers::checkGroupConstraintsProd( - group.name(), - group.parent(), - parentGroup, - this->wellState(), - this->groupState(), - reportStepIdx, - guideRate_.get(), - rates.data(), - phase_usage_, - group.getGroupEfficiencyFactor(), - schedule(), - summaryState, - resv_coeff, - deferred_logger); - if (changed.first) { - switched_groups.insert(group.name()); - const auto exceed_action = group.productionControls(summaryState).exceed_action; - actionOnBrokenConstraints(group, exceed_action, Group::ProductionCMode::FLD, deferred_logger); - } - } - } - - // call recursively down the group hiearchy - for (const std::string& groupName : group.groups()) { - checkGroupHigherConstraints( schedule().getGroup(groupName, reportStepIdx), deferred_logger, switched_groups); - } - } - - - - - - - - - template - void - BlackoilWellModel:: - updateEclWells(const int timeStepIdx, const std::unordered_set& wells) { - const auto& schedule = this->ebosSimulator_.vanguard().schedule(); - for (const auto& wname : wells) { - 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()) { - auto well_index = std::distance( this->wells_ecl_.begin(), well_iter ); - this->wells_ecl_[well_index] = schedule.getWell(wname, timeStepIdx); - - const auto& well = this->wells_ecl_[well_index]; - auto& pd = this->well_perf_data_[well_index]; - auto pdIter = pd.begin(); - for (const auto& conn : well.getConnections()) { - if (conn.state() != Connection::State::SHUT) { - pdIter->connection_transmissibility_factor = conn.CF(); - ++pdIter; - } - } - this->wellState().updateStatus(well_index, well.getStatus()); - this->wellState().resetConnectionTransFactors(well_index, pd); - this->prod_index_calc_[well_index].reInit(well); - } - } - } - - - - template - double - BlackoilWellModel:: - wellPI(const int well_index) const - { - const auto& pu = this->phase_usage_; - const auto& pi = this->wellState().productivityIndex(well_index); - - const auto preferred = this->wells_ecl_[well_index].getPreferredPhase(); - switch (preferred) { // Should really have LIQUID = OIL + WATER here too... - case Phase::WATER: - return pu.phase_used[BlackoilPhases::PhaseIndex::Aqua] - ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Aqua]] - : 0.0; - - case Phase::OIL: - return pu.phase_used[BlackoilPhases::PhaseIndex::Liquid] - ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Liquid]] - : 0.0; - - case Phase::GAS: - return pu.phase_used[BlackoilPhases::PhaseIndex::Vapour] - ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Vapour]] - : 0.0; - - default: - throw std::invalid_argument { - "Unsupported preferred phase " + - 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_index = std::distance(this->wells_ecl_.begin(), well_iter); - return this->wellPI(well_index); - } - - - - - template int BlackoilWellModel:: @@ -2853,6 +1712,27 @@ namespace Opm { + template + void + BlackoilWellModel:: + setWellWsolvent(const std::string& wellName, double wsolvent) { + auto well = getWell(wellName); + well->setWsolvent(wsolvent); + } + + + + template + void + BlackoilWellModel:: + calcRates(const int fipnum, + const int pvtreg, + std::vector& resv_coeff) + { + rateConverter_->calcCoeff(fipnum, pvtreg, resv_coeff); + } + + template void BlackoilWellModel:: @@ -2935,414 +1815,6 @@ 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:: - updateWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, const WellState& wellState) { - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateWsolvent(groupTmp, schedule, reportStepIdx, wellState); - } - - if (group.isProductionGroup()) - return; - - auto currentGroupControl = this->groupState().injection_control(group.name(), Phase::GAS); - if( currentGroupControl == Group::InjectionCMode::REIN ) { - int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; - const auto& summaryState = ebosSimulator_.vanguard().summaryState(); - const auto& controls = group.injectionControls(Phase::GAS, summaryState); - const Group& groupRein = schedule.getGroup(controls.reinj_group, reportStepIdx); - double gasProductionRate = WellGroupHelpers::sumWellRates(groupRein, schedule, wellState, reportStepIdx, gasPos, /*isInjector*/false); - double solventProductionRate = WellGroupHelpers::sumSolventRates(groupRein, schedule, wellState, reportStepIdx, /*isInjector*/false); - - const auto& comm = ebosSimulator_.vanguard().grid().comm(); - solventProductionRate = comm.sum(solventProductionRate); - gasProductionRate = comm.sum(gasProductionRate); - - double wsolvent = 0.0; - if (std::abs(gasProductionRate) > 1e-6) - wsolvent = solventProductionRate / gasProductionRate; - - setWsolvent(group, schedule, reportStepIdx, wsolvent); - } - } - - template - void - BlackoilWellModel:: - setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent) { - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - setWsolvent(groupTmp, schedule, reportStepIdx, wsolvent); - } - - for (const std::string& wellName : group.wells()) { - const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); - if (wellTmp.getStatus() == Well::Status::SHUT) - continue; - - auto well = getWell(wellName); - well->setWsolvent(wsolvent); - } - } - - - - - - template - void - BlackoilWellModel:: - assignWellGuideRates(data::Wells& wsrpt) const - { - for (const auto& well : this->wells_ecl_) { - auto xwPos = wsrpt.find(well.name()); - if (xwPos == wsrpt.end()) { // No well results. Unexpected. - continue; - } - - xwPos->second.guide_rates = this->getGuideRateValues(well); - } - } - - - - - - template - void - BlackoilWellModel:: - assignShutConnections(data::Wells& wsrpt) const - { - auto wellID = 0; - - for (const auto& well : this->wells_ecl_) { - 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 = xwel.connections; - for (const auto& conn : well.getConnections()) { - if (skip(conn)) { - continue; - } - - auto& xc = xcon.emplace_back(); - xc.index = conn.global_index(); - xc.pressure = xc.reservoir_rate = 0.0; - - xc.effective_Kh = conn.Kh(); - xc.trans_factor = conn.CF(); - } - - ++wellID; - } - } - - - - - - template - void - BlackoilWellModel:: - assignGroupValues(const int reportStepIdx, - const Schedule& sched, - std::map& gvalues) const - { - const auto groupGuideRates = - this->calculateAllGroupGuiderates(reportStepIdx, sched); - - for (const auto& gname : sched.groupNames(reportStepIdx)) { - const auto& grup = sched.getGroup(gname, reportStepIdx); - - auto& gdata = gvalues[gname]; - this->assignGroupControl(grup, gdata); - this->assignGroupGuideRates(grup, groupGuideRates, gdata); - } - } - - template - void - BlackoilWellModel:: - assignNodeValues(std::map& nodevalues) const - { - nodevalues.clear(); - for (const auto& [node, pressure] : node_pressures_) { - nodevalues.emplace(node, data::NodeData{pressure}); - } - } - - template - std::unordered_map - BlackoilWellModel:: - calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const - { - auto gr = std::unordered_map{}; - auto up = std::vector{}; - - // Start at well level, accumulate contributions towards root of - // group tree (FIELD group). - - for (const auto& wname : sched.wellNames(reportStepIdx)) { - if (! (this->wellState().hasWellRates(wname) && - this->guideRate_->has(wname))) - { - continue; - } - - const auto& well = sched.getWell(wname, reportStepIdx); - const auto& parent = well.groupName(); - - if (parent == "FIELD") { - // Well parented directly to "FIELD". Inadvisable and - // unexpected, but nothing to do about that here. Just skip - // this guide rate contribution. - continue; - } - - auto& grval = well.isInjector() - ? gr[parent].injection - : gr[parent].production; - - grval += this->getGuideRateValues(well); - up.push_back(parent); - } - - // Propagate accumulated guide rates up towards root of group tree. - // Override accumulation if there is a GUIDERAT specification that - // applies to a group. - std::sort(up.begin(), up.end()); - auto start = 0*up.size(); - auto u = std::unique(up.begin(), up.end()); - auto nu = std::distance(up.begin(), u); - while (nu > 0) { - const auto ntot = up.size(); - - for (auto gi = 0*nu; gi < nu; ++gi) { - const auto& gname = up[start + gi]; - const auto& group = sched.getGroup(gname, reportStepIdx); - - if (this->guideRate_->has(gname)) { - gr[gname].production = this->getGuideRateValues(group); - } - - if (this->guideRate_->has(gname, Phase::WATER) - || this->guideRate_->has(gname, Phase::GAS)) { - gr[gname].injection = this->getGuideRateInjectionGroupValues(group); - } - - const auto parent = group.parent(); - if (parent == "FIELD") { continue; } - - gr[parent].injection += gr[gname].injection; - gr[parent].production += gr[gname].production; - up.push_back(parent); - } - - start = ntot; - - auto begin = up.begin() + ntot; - std::sort(begin, up.end()); - u = std::unique(begin, up.end()); - nu = std::distance(begin, u); - } - - return gr; - } - - template - void - BlackoilWellModel:: - assignGroupControl(const Group& group, data::GroupData& gdata) const - { - const auto& gname = group.name(); - const auto grup_type = group.getGroupType(); - auto& cgc = gdata.currentControl; - - cgc.currentProdConstraint = - ::Opm::Group::ProductionCMode::NONE; - - cgc.currentGasInjectionConstraint = - cgc.currentWaterInjectionConstraint = - ::Opm::Group::InjectionCMode::NONE; - - if (this->groupState().has_production_control(gname)) { - cgc.currentProdConstraint = this->groupState().production_control(gname); - } - - if ((grup_type == ::Opm::Group::GroupType::INJECTION) || - (grup_type == ::Opm::Group::GroupType::MIXED)) - { - if (this->groupState().has_injection_control(gname, ::Opm::Phase::WATER)) { - cgc.currentWaterInjectionConstraint = this->groupState().injection_control(gname, Phase::WATER); - } - - if (this->groupState().has_injection_control(gname, ::Opm::Phase::GAS)) { - cgc.currentGasInjectionConstraint = this->groupState().injection_control(gname, Phase::GAS); - } - } - } - - template - data::GuideRateValue - BlackoilWellModel:: - getGuideRateValues(const Well& well) const - { - auto grval = data::GuideRateValue{}; - - assert (this->guideRate_ != nullptr); - - const auto& wname = well.name(); - if (! this->wellState().hasWellRates(wname)) { - // No flow rates for 'wname' -- might be before well comes - // online (e.g., for the initial condition before simulation - // starts). - return grval; - } - - if (! this->guideRate_->has(wname)) { - // No guiderates exist for 'wname'. - return grval; - } - - const auto qs = WellGroupHelpers:: - getWellRateVector(this->wellState(), this->phase_usage_, wname); - - this->getGuideRateValues(qs, well.isInjector(), wname, grval); - - return grval; - } - template - data::GuideRateValue - BlackoilWellModel:: - getGuideRateInjectionGroupValues(const Group& group) const - { - auto grval = data::GuideRateValue{}; - - assert (this->guideRate_ != nullptr); - - const auto& gname = group.name(); - if (this->guideRate_->has(gname, Phase::GAS)) { - grval.set(data::GuideRateValue::Item::Gas, - this->guideRate_->get(gname, Phase::GAS)); - } - if (this->guideRate_->has(gname, Phase::WATER)) { - grval.set(data::GuideRateValue::Item::Water, - this->guideRate_->get(gname, Phase::WATER)); - } - return grval; - } - - template - data::GuideRateValue - BlackoilWellModel:: - getGuideRateValues(const Group& group) const - { - auto grval = data::GuideRateValue{}; - - assert (this->guideRate_ != nullptr); - - const auto& gname = group.name(); - - if ( ! this->groupState().has_production_rates(gname)) { - // No flow rates for production group 'gname' -- might be before group comes - // online (e.g., for the initial condition before simulation - // starts). - return grval; - } - - if (! this->guideRate_->has(gname)) { - // No guiderates exist for 'gname'. - return grval; - } - - const auto qs = WellGroupHelpers::getProductionGroupRateVector(this->groupState(), this->phase_usage_, gname); - - const auto is_inj = false; // This procedure only applies to G*PGR. - this->getGuideRateValues(qs, is_inj, gname, grval); - - return grval; - } - - template - void - BlackoilWellModel:: - getGuideRateValues(const GuideRate::RateVector& qs, - const bool is_inj, - const std::string& wgname, - data::GuideRateValue& grval) const - { - auto getGR = [this, &wgname, &qs](const GuideRateModel::Target t) - { - return this->guideRate_->get(wgname, t, qs); - }; - - // Note: GuideRate does currently (2020-07-20) not support Target::RES. - grval.set(data::GuideRateValue::Item::Gas, - getGR(GuideRateModel::Target::GAS)); - - grval.set(data::GuideRateValue::Item::Water, - getGR(GuideRateModel::Target::WAT)); - - if (! is_inj) { - // Producer. Extract "all" guiderate values. - grval.set(data::GuideRateValue::Item::Oil, - getGR(GuideRateModel::Target::OIL)); - } - } - - template - void - BlackoilWellModel:: - assignGroupGuideRates(const Group& group, - const std::unordered_map& groupGuideRates, - data::GroupData& gdata) const - { - auto& prod = gdata.guideRates.production; prod.clear(); - auto& inj = gdata.guideRates.injection; inj .clear(); - - auto xgrPos = groupGuideRates.find(group.name()); - if ((xgrPos == groupGuideRates.end()) || - !this->guideRate_->has(group.name())) - { - // No guiderates defined for this group. - return; - } - - const auto& xgr = xgrPos->second; - prod = xgr.production; - inj = xgr.injection; - } - template void BlackoilWellModel:: diff --git a/tests/test_glift1.cpp b/tests/test_glift1.cpp index 2425dc037..d4235b472 100644 --- a/tests/test_glift1.cpp +++ b/tests/test_glift1.cpp @@ -160,7 +160,7 @@ BOOST_AUTO_TEST_CASE(G1) const auto &well = wells_ecl[*idx]; BOOST_CHECK_EQUAL( well.name(), "B-1H"); const auto& summary_state = simulator->vanguard().summaryState(); - WellState &well_state = const_cast(well_model.wellState()); + WellState &well_state = well_model.wellState(); GasLiftSingleWell glift {*std_well, *(simulator.get()), summary_state, deferred_logger, well_state}; auto state = glift.runOptimize(simulator->model().newtonMethod().numIterations());