23#ifndef OPM_BLACKOILWELLMODEL_IMPL_HEADER_INCLUDED
24#define OPM_BLACKOILWELLMODEL_IMPL_HEADER_INCLUDED
27#ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
32#include <opm/grid/utility/cartesianToCompressed.hpp>
33#include <opm/common/utility/numeric/RootFinders.hpp>
35#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
36#include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
37#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
38#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
39#include <opm/input/eclipse/Schedule/Well/WellTestConfig.hpp>
41#include <opm/input/eclipse/Units/UnitSystem.hpp>
68#include <fmt/format.h>
71 template<
typename TypeTag>
74 : WellConnectionModule(*this, simulator.gridView().comm())
75 , BlackoilWellModelGeneric<Scalar>(simulator.vanguard().schedule(),
77 simulator.vanguard().summaryState(),
78 simulator.vanguard().eclState(),
80 simulator.gridView().comm())
81 , simulator_(simulator)
82 , guide_rate_handler_{
84 simulator.vanguard().schedule(),
85 simulator.vanguard().summaryState(),
86 simulator.vanguard().grid().comm()
88 , gaslift_(this->terminal_output_, this->phase_usage_)
90 local_num_cells_ = simulator_.gridView().size(0);
93 global_num_cells_ = simulator_.vanguard().globalNumCells();
96 auto& parallel_wells = simulator.vanguard().parallelWells();
98 this->parallel_well_info_.reserve(parallel_wells.size());
99 for(
const auto& name_bool : parallel_wells) {
100 this->parallel_well_info_.emplace_back(name_bool, grid().comm());
104 this->alternative_well_rate_init_ =
105 Parameters::Get<Parameters::AlternativeWellRateInit>();
107 using SourceDataSpan =
108 typename PAvgDynamicSourceData<Scalar>::template SourceDataSpan<Scalar>;
110 this->wbp_.initializeSources(
111 [
this](
const std::size_t globalIndex)
112 {
return this->compressedIndexForInterior(globalIndex); },
113 [
this](
const int localCell, SourceDataSpan sourceTerms)
115 using Item =
typename SourceDataSpan::Item;
117 const auto* intQuants = this->simulator_.model()
118 .cachedIntensiveQuantities(localCell, 0);
119 const auto& fs = intQuants->fluidState();
122 .set(Item::PoreVol, intQuants->porosity().value() *
123 this->simulator_.model().dofTotalVolume(localCell))
124 .set(Item::Depth, this->depth_[localCell]);
126 constexpr auto io = FluidSystem::oilPhaseIdx;
127 constexpr auto ig = FluidSystem::gasPhaseIdx;
128 constexpr auto iw = FluidSystem::waterPhaseIdx;
131 const auto haveOil = FluidSystem::phaseIsActive(io);
132 const auto haveGas = FluidSystem::phaseIsActive(ig);
133 const auto haveWat = FluidSystem::phaseIsActive(iw);
135 auto weightedPhaseDensity = [&fs](
const auto ip)
137 return fs.saturation(ip).value() * fs.density(ip).value();
140 if (haveOil) { sourceTerms.set(Item::Pressure, fs.pressure(io).value()); }
141 else if (haveGas) { sourceTerms.set(Item::Pressure, fs.pressure(ig).value()); }
142 else { sourceTerms.set(Item::Pressure, fs.pressure(iw).value()); }
146 if (haveOil) { rho += weightedPhaseDensity(io); }
147 if (haveGas) { rho += weightedPhaseDensity(ig); }
148 if (haveWat) { rho += weightedPhaseDensity(iw); }
150 sourceTerms.set(Item::MixtureDensity, rho);
155 template<
typename TypeTag>
162 template<
typename TypeTag>
167 extractLegacyCellPvtRegionIndex_();
168 extractLegacyDepth_();
170 gravity_ = simulator_.problem().gravity()[2];
172 this->initial_step_ =
true;
175 simulator_.model().addAuxiliaryModule(
this);
177 is_cell_perforated_.resize(local_num_cells_,
false);
181 template<
typename TypeTag>
186 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
187 + ScheduleEvents::NEW_WELL;
188 const auto& events = this->schedule()[reportStepIdx].wellgroup_events();
189 for (
auto& wellPtr : this->well_container_) {
190 const bool well_opened_this_step = this->report_step_starts_ &&
191 events.hasEvent(wellPtr->name(),
192 effective_events_mask);
193 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
194 this->B_avg_, well_opened_this_step);
198 template<
typename TypeTag>
205 this->report_step_starts_ =
true;
206 this->report_step_start_events_ = this->schedule()[timeStepIdx].wellgroup_events();
208 this->rateConverter_ = std::make_unique<RateConverterType>
209 (this->phase_usage_, std::vector<int>(this->local_num_cells_, 0));
213 const auto enableWellPIScaling =
true;
214 this->initializeLocalWellStructure(timeStepIdx, enableWellPIScaling);
217 this->initializeGroupStructure(timeStepIdx);
219 const auto& comm = this->simulator_.vanguard().grid().comm();
225 this->rateConverter_->template defineState<ElementContext>(this->simulator_);
229 const auto& sched_state = this->schedule()[timeStepIdx];
231 this->vfp_properties_ = std::make_unique<VFPProperties<Scalar>>
232 (sched_state.vfpinj(), sched_state.vfpprod(), this->wellState());
236 "beginReportStep() failed: ",
237 this->terminal_output_, comm)
241 this->commitWGState();
243 this->wellStructureChangedDynamically_ =
false;
250 template <
typename TypeTag>
254 const bool enableWellPIScaling)
258 const auto& comm = this->simulator_.vanguard().grid().comm();
261 this->wells_ecl_ = this->getLocalWells(reportStepIdx);
262 this->local_parallel_well_info_ =
263 this->createLocalParallelWellInfo(this->wells_ecl_);
270 this->initializeWellPerfData();
271 this->initializeWellState(reportStepIdx);
272 this->wbp_.initializeWBPCalculationService();
274 if (this->param_.use_multisegment_well_ && this->anyMSWellOpenLocal()) {
275 this->wellState().initWellStateMSWell(this->wells_ecl_, &this->prevWellState());
278 this->initializeWellProdIndCalculators();
280 if (enableWellPIScaling && this->schedule()[reportStepIdx].events()
281 .hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX))
283 this->runWellPIScaling(reportStepIdx, local_deferredLogger);
287 "Failed to initialize local well structure: ",
288 this->terminal_output_, comm)
295 template <
typename TypeTag>
302 const auto& comm = this->simulator_.vanguard().grid().comm();
306 const auto& fieldGroup =
307 this->schedule().getGroup(
"FIELD", reportStepIdx);
311 this->summaryState(),
317 if (this->schedule()[reportStepIdx].has_gpmaint()) {
322 this->eclState_.fieldProps(),
324 this->regionalAveragePressureCalculator_);
328 "Failed to initialize group structure: ",
329 this->terminal_output_, comm)
337 template<
typename TypeTag>
342 OPM_TIMEBLOCK(beginTimeStep);
344 this->updateAverageFormationFactor();
348 this->switched_prod_groups_.clear();
349 this->switched_inj_groups_.clear();
351 if (this->wellStructureChangedDynamically_) {
356 const auto reportStepIdx =
357 this->simulator_.episodeIndex();
361 const auto enableWellPIScaling =
false;
363 this->initializeLocalWellStructure(reportStepIdx, enableWellPIScaling);
364 this->initializeGroupStructure(reportStepIdx);
366 this->commitWGState();
372 this->wellStructureChangedDynamically_ =
false;
375 this->resetWGState();
377 const int reportStepIdx = simulator_.episodeIndex();
378 this->updateAndCommunicateGroupData(reportStepIdx,
379 simulator_.model().newtonMethod().numIterations(),
380 param_.nupcol_group_rate_tolerance_,
false,
381 local_deferredLogger);
383 this->wellState().updateWellsDefaultALQ(this->schedule(), reportStepIdx, this->summaryState());
384 this->wellState().gliftTimeStepInit();
386 const double simulationTime = simulator_.time();
390 wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
393 createWellContainer(reportStepIdx);
396 const Grid& grid = simulator_.vanguard().grid();
397 this->wells_active_ = grid.comm().max(!this->well_container_.empty());
402 this->initWellContainer(reportStepIdx);
405 std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(),
false);
406 for (
auto& well : well_container_) {
407 well->updatePerforatedCell(is_cell_perforated_);
411 this->calculateEfficiencyFactors(reportStepIdx);
413 if constexpr (has_polymer_)
415 if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
416 this->setRepRadiusPerfLength();
423 this->terminal_output_, simulator_.vanguard().grid().comm());
425 for (
auto& well : well_container_) {
426 well->setVFPProperties(this->vfp_properties_.get());
427 well->setGuideRate(&this->guideRate_);
430 this->updateFiltrationModelsPreStep(local_deferredLogger);
433 for (
auto& well : well_container_) {
434 well->closeCompletions(this->wellTestState());
440 if (alternative_well_rate_init_) {
445 for (
const auto& well : well_container_) {
446 if (well->isProducer() && !well->wellIsStopped()) {
447 well->initializeProducerWellState(simulator_, this->wellState(), local_deferredLogger);
452 for (
const auto& well : well_container_) {
453 if (well->isVFPActive(local_deferredLogger)){
454 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
458 this->updateWellPotentials(reportStepIdx,
460 simulator_.vanguard().summaryConfig(),
461 local_deferredLogger);
462 }
catch ( std::runtime_error& e ) {
463 const std::string msg =
"A zero well potential is returned for output purposes. ";
464 local_deferredLogger.
warning(
"WELL_POTENTIAL_CALCULATION_FAILED", msg);
466 this->guide_rate_handler_.setLogger(&local_deferredLogger);
467#ifdef RESERVOIR_COUPLING_ENABLED
468 if (this->isReservoirCouplingMaster()) {
469 this->guide_rate_handler_.receiveMasterGroupPotentialsFromSlaves();
473 this->guide_rate_handler_.updateGuideRates(
474 reportStepIdx, simulationTime, this->wellState(), this->groupState()
476#ifdef RESERVOIR_COUPLING_ENABLED
477 if (this->isReservoirCouplingSlave()) {
478 this->guide_rate_handler_.sendSlaveGroupPotentialsToMaster(this->groupState());
484 if (this->schedule_[reportStepIdx].has_gpmaint()) {
485 for (
const auto& calculator : regionalAveragePressureCalculator_) {
486 calculator.second->template defineState<ElementContext>(simulator_);
488 const double dt = simulator_.timeStepSize();
489 const Group& fieldGroup = this->schedule().getGroup(
"FIELD", reportStepIdx);
492 regionalAveragePressureCalculator_,
499 this->updateAndCommunicateGroupData(reportStepIdx,
500 simulator_.model().newtonMethod().numIterations(),
501 param_.nupcol_group_rate_tolerance_,
503 local_deferredLogger);
506 for (
auto& well : well_container_) {
507 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
508 + ScheduleEvents::INJECTION_TYPE_CHANGED
509 + ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
510 + ScheduleEvents::NEW_WELL;
512 const auto& events = this->schedule()[reportStepIdx].wellgroup_events();
513 const bool event = this->report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
514 const bool dyn_status_change = this->wellState().well(well->name()).status
515 != this->prevWellState().well(well->name()).status;
517 if (event || dyn_status_change) {
519 well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), local_deferredLogger);
520 well->calculateExplicitQuantities(simulator_, this->wellState(), local_deferredLogger);
521 well->solveWellEquation(simulator_, this->wellState(), this->groupState(), local_deferredLogger);
522 }
catch (
const std::exception& e) {
523 const std::string msg =
"Compute initial well solution for new well " + well->name() +
" failed. Continue with zero initial rates";
524 local_deferredLogger.
warning(
"WELL_INITIAL_SOLVE_FAILED", msg);
533 const std::string msg =
"Compute initial well solution for new wells failed. Continue with zero initial rates";
534 local_deferredLogger.
warning(
"WELL_INITIAL_SOLVE_FAILED", msg);
537 const auto& comm = simulator_.vanguard().grid().comm();
539 exc_type,
"beginTimeStep() failed: " + exc_msg, this->terminal_output_, comm);
543 template<
typename TypeTag>
546 const double simulationTime,
549 for (
const std::string& well_name : this->getWellsForTesting(timeStepIdx, simulationTime)) {
550 const Well& wellEcl = this->schedule().getWell(well_name, timeStepIdx);
551 if (wellEcl.getStatus() == Well::Status::SHUT)
554 WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
556 well->init(&this->phase_usage_, depth_, gravity_, B_avg_,
true);
558 Scalar well_efficiency_factor = wellEcl.getEfficiencyFactor() *
559 this->wellState().getGlobalEfficiencyScalingFactor(well_name);
564 well_efficiency_factor);
566 well->setWellEfficiencyFactor(well_efficiency_factor);
567 well->setVFPProperties(this->vfp_properties_.get());
568 well->setGuideRate(&this->guideRate_);
571 if (well->isProducer()) {
572 well->initializeProducerWellState(simulator_, this->wellState(), deferred_logger);
574 if (well->isVFPActive(deferred_logger)) {
575 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
578 const auto& network = this->schedule()[timeStepIdx].network();
579 if (network.active() && !this->node_pressures_.empty()) {
580 if (well->isProducer()) {
581 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
582 if (it != this->node_pressures_.end()) {
585 const Scalar nodal_pressure = it->second;
586 well->setDynamicThpLimit(nodal_pressure);
592 GLiftEclWells ecl_well_map;
593 gaslift_.initGliftEclWellMap(well_container_, ecl_well_map);
594 well->wellTesting(simulator_,
598 this->wellTestState(),
601 this->well_open_times_,
603 }
catch (
const std::exception& e) {
604 const std::string msg = fmt::format(
"Exception during testing of well: {}. The well will not open.\n Exception message: {}", wellEcl.name(), e.what());
605 deferred_logger.
warning(
"WELL_TESTING_FAILED", msg);
615 template<
typename TypeTag>
621 for (
auto&& pinfo : this->local_parallel_well_info_)
632 template<
typename TypeTag>
642 template<
typename TypeTag>
647 this->closed_this_step_.clear();
650 this->report_step_starts_ =
false;
651 const int reportStepIdx = simulator_.episodeIndex();
654 for (
const auto& well : well_container_) {
655 if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
656 well->updateWaterThroughput(dt, this->wellState());
660 for (
const auto& well : well_container_) {
661 well->updateConnectionTransmissibilityFactor(simulator_, this->wellState().well(well->indexOfWell()));
662 well->updateConnectionDFactor(simulator_, this->wellState().well(well->indexOfWell()));
665 if (Indices::waterEnabled) {
666 this->updateFiltrationModelsPostStep(dt, FluidSystem::waterPhaseIdx, local_deferredLogger);
670 this->updateInjMult(local_deferredLogger);
673 for (
const auto& well : well_container_) {
674 well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);
677 if (this->terminal_output_) {
678 this->reportGroupSwitching(local_deferredLogger);
682 rateConverter_->template defineState<ElementContext>(simulator_);
686 this->updateWellPotentials(reportStepIdx,
688 simulator_.vanguard().summaryConfig(),
689 local_deferredLogger);
690 }
catch ( std::runtime_error& e ) {
691 const std::string msg =
"A zero well potential is returned for output purposes. ";
692 local_deferredLogger.
warning(
"WELL_POTENTIAL_CALCULATION_FAILED", msg);
695 updateWellTestState(simulationTime, this->wellTestState());
698 const Group& fieldGroup = this->schedule_.getGroup(
"FIELD", reportStepIdx);
699 this->checkGEconLimits(fieldGroup, simulationTime,
700 simulator_.episodeIndex(), local_deferredLogger);
701 this->checkGconsaleLimits(fieldGroup, this->wellState(),
702 simulator_.episodeIndex(), local_deferredLogger);
704 this->calculateProductivityIndexValues(local_deferredLogger);
706 this->commitWGState();
710 if (this->terminal_output_) {
715 this->computeWellTemperature();
719 template<
typename TypeTag>
723 unsigned elemIdx)
const
727 if (!is_cell_perforated_[elemIdx]) {
731 for (
const auto& well : well_container_)
732 well->addCellRates(rate, elemIdx);
736 template<
typename TypeTag>
737 template <
class Context>
741 const Context& context,
743 unsigned timeIdx)
const
746 int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
748 if (!is_cell_perforated_[elemIdx]) {
752 for (
const auto& well : well_container_)
753 well->addCellRates(rate, elemIdx);
758 template<
typename TypeTag>
763 const auto pressIx = []()
765 if (Indices::oilEnabled) {
return FluidSystem::oilPhaseIdx; }
766 if (Indices::waterEnabled) {
return FluidSystem::waterPhaseIdx; }
768 return FluidSystem::gasPhaseIdx;
771 auto cellPressures = std::vector<Scalar>(this->local_num_cells_,
Scalar{0});
772 auto cellTemperatures = std::vector<Scalar>(this->local_num_cells_,
Scalar{0});
775 const auto& gridView = this->simulator_.vanguard().gridView();
778 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
779 elemCtx.updatePrimaryStencil(elem);
780 elemCtx.updatePrimaryIntensiveQuantities(0);
782 const auto ix = elemCtx.globalSpaceIndex(0, 0);
783 const auto& fs = elemCtx.intensiveQuantities(0, 0).fluidState();
785 cellPressures[ix] = fs.pressure(pressIx).value();
786 cellTemperatures[ix] = fs.temperature(0).value();
789 this->simulator_.vanguard().grid().comm());
791 this->wellState().init(cellPressures, cellTemperatures, this->schedule(), this->wells_ecl_,
792 this->local_parallel_well_info_, timeStepIdx,
793 &this->prevWellState(), this->well_perf_data_,
794 this->summaryState(), simulator_.vanguard().enableDistributedWells());
801 template<
typename TypeTag>
808 const int nw = this->numLocalWells();
810 well_container_.clear();
813 well_container_.reserve(nw);
815 const auto& wmatcher = this->schedule().wellMatcher(report_step);
816 const auto& wcycle = this->schedule()[report_step].wcycle.get();
820 std::for_each(this->wells_ecl_.begin(), this->wells_ecl_.end(),
821 [
this, &wg_events = this->report_step_start_events_](
const auto& well_ecl)
823 if (!well_ecl.hasConnections()) {
828 constexpr auto events_mask = ScheduleEvents::WELL_STATUS_CHANGE |
829 ScheduleEvents::REQUEST_OPEN_WELL;
830 const bool well_status_change =
831 this->report_step_starts_ &&
832 wg_events.hasEvent(well_ecl.name(), events_mask);
833 if (well_status_change) {
834 if (well_ecl.getStatus() == WellStatus::OPEN) {
835 this->well_open_times_.insert_or_assign(well_ecl.name(),
836 this->simulator_.time());
837 this->well_close_times_.erase(well_ecl.name());
838 }
else if (well_ecl.getStatus() == WellStatus::SHUT) {
839 this->well_close_times_.insert_or_assign(well_ecl.name(),
840 this->simulator_.time());
841 this->well_open_times_.erase(well_ecl.name());
847 const auto cycle_states = wcycle.wellStatus(this->simulator_.time(),
849 this->well_open_times_,
850 this->well_close_times_);
852 for (
int w = 0; w < nw; ++w) {
853 const Well& well_ecl = this->wells_ecl_[w];
855 if (!well_ecl.hasConnections()) {
860 const std::string& well_name = well_ecl.name();
861 const auto well_status = this->schedule()
862 .getWell(well_name, report_step).getStatus();
864 if ((well_ecl.getStatus() == Well::Status::SHUT) ||
865 (well_status == Well::Status::SHUT))
868 if (well_ecl.getStatus() != Well::Status::SHUT) {
869 this->closed_this_step_.insert(well_name);
870 this->wellState().shutWell(w);
873 this->well_open_times_.erase(well_name);
874 this->well_close_times_.erase(well_name);
879 if (this->wellTestState().well_is_closed(well_name)) {
884 const bool closed_this_step = (this->wellTestState().lastTestTime(well_name) == simulator_.time());
887 auto& events = this->wellState().well(w).events;
888 if (events.hasEvent(ScheduleEvents::REQUEST_OPEN_WELL)) {
889 if (!closed_this_step) {
890 this->wellTestState().open_well(well_name);
891 this->wellTestState().open_completions(well_name);
892 this->well_open_times_.insert_or_assign(well_name,
893 this->simulator_.time());
894 this->well_close_times_.erase(well_name);
896 events.clearEvent(ScheduleEvents::REQUEST_OPEN_WELL);
902 bool wellIsStopped =
false;
903 if (this->wellTestState().well_is_closed(well_name))
905 if (well_ecl.getAutomaticShutIn()) {
907 this->wellState().shutWell(w);
908 this->well_close_times_.erase(well_name);
909 this->well_open_times_.erase(well_name);
912 if (!well_ecl.getAllowCrossFlow()) {
915 this->wellState().shutWell(w);
916 this->well_close_times_.erase(well_name);
917 this->well_open_times_.erase(well_name);
921 this->wellState().stopWell(w);
922 wellIsStopped =
true;
927 if (!well_ecl.getAllowCrossFlow()) {
928 const bool any_zero_rate_constraint = well_ecl.isProducer()
929 ? well_ecl.productionControls(this->summaryState_).anyZeroRateConstraint()
930 : well_ecl.injectionControls(this->summaryState_).anyZeroRateConstraint();
931 if (any_zero_rate_constraint) {
933 local_deferredLogger.
debug(fmt::format(
" Well {} gets shut due to having zero rate constraint and disallowing crossflow ", well_ecl.name()) );
934 this->wellState().shutWell(w);
935 this->well_close_times_.erase(well_name);
936 this->well_open_times_.erase(well_name);
941 if (well_status == Well::Status::STOP) {
942 this->wellState().stopWell(w);
943 this->well_close_times_.erase(well_name);
944 this->well_open_times_.erase(well_name);
945 wellIsStopped =
true;
948 if (!wcycle.empty()) {
949 const auto it = cycle_states.find(well_name);
950 if (it != cycle_states.end()) {
952 this->wellState().shutWell(w);
955 this->wellState().openWell(w);
960 well_container_.emplace_back(this->createWellPointer(w, report_step));
963 well_container_.back()->stopWell();
964 this->well_close_times_.erase(well_name);
965 this->well_open_times_.erase(well_name);
969 if (!wcycle.empty()) {
970 const auto schedule_open =
971 [&wg_events = this->report_step_start_events_](
const std::string& name)
973 return wg_events.hasEvent(name, ScheduleEvents::REQUEST_OPEN_WELL);
975 for (
const auto& [wname, wscale] : wcycle.efficiencyScale(this->simulator_.time(),
976 this->simulator_.timeStepSize(),
978 this->well_open_times_,
981 this->wellState().updateEfficiencyScalingFactor(wname, wscale);
982 this->schedule_.add_event(ScheduleEvents::WELLGROUP_EFFICIENCY_UPDATE, report_step);
991 if (this->terminal_output_) {
995 this->well_container_generic_.clear();
996 for (
auto& w : well_container_)
997 this->well_container_generic_.push_back(w.get());
999 const auto& network = this->schedule()[report_step].network();
1000 if (network.active() && !this->node_pressures_.empty()) {
1001 for (
auto& well: this->well_container_generic_) {
1005 if (well->isProducer()) {
1006 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
1007 if (it != this->node_pressures_.end()) {
1010 const Scalar nodal_pressure = it->second;
1011 well->setDynamicThpLimit(nodal_pressure);
1017 this->wbp_.registerOpenWellsForWBPCalculation();
1024 template <
typename TypeTag>
1025 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1029 const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
1031 if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
1032 return this->
template createTypedWellPointer<StandardWell<TypeTag>>(wellID, report_step);
1035 return this->
template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, report_step);
1043 template <
typename TypeTag>
1044 template <
typename WellType>
1045 std::unique_ptr<WellType>
1050 const auto& perf_data = this->well_perf_data_[wellID];
1053 const auto pvtreg = perf_data.
empty()
1054 ? 0 : this->pvt_region_idx_[perf_data.front().cell_index];
1056 const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get();
1057 const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
1059 return std::make_unique<WellType>(this->wells_ecl_[wellID],
1063 *this->rateConverter_,
1065 this->numComponents(),
1075 template<
typename TypeTag>
1079 const int report_step,
1083 const auto it = std::find_if(this->wells_ecl_.begin(),
1084 this->wells_ecl_.end(),
1085 [&well_name](
const auto& w)
1086 { return well_name == w.name(); });
1088 if (it == this->wells_ecl_.end()) {
1090 fmt::format(
"Could not find well {} in wells_ecl ", well_name),
1094 const int pos =
static_cast<int>(std::distance(this->wells_ecl_.begin(), it));
1095 return this->createWellPointer(pos, report_step);
1100 template<
typename TypeTag>
1106 const double dt = this->simulator_.timeStepSize();
1108 auto& well_state = this->wellState();
1110 const bool changed_well_group = updateWellControlsAndNetwork(
true, dt, deferred_logger);
1111 assembleWellEqWithoutIteration(dt, deferred_logger);
1112 const bool converged = this->getWellConvergence(this->B_avg_,
true).converged() && !changed_well_group;
1115 for (
auto& well : this->well_container_) {
1116 well->solveEqAndUpdateWellState(simulator_, well_state, deferred_logger);
1119 this->simulator_.vanguard().grid().comm());
1122 const std::string msg = fmt::format(
"Initial (pre-step) network balance did not converge.");
1130 template<
typename TypeTag>
1139 this->guide_rate_handler_.setLogger(&local_deferredLogger);
1141 if (gaslift_.terminalOutput()) {
1142 const std::string msg =
1143 fmt::format(
"assemble() : iteration {}" , iterationIdx);
1144 gaslift_.gliftDebug(msg, local_deferredLogger);
1148 Dune::Timer perfTimer;
1150 this->closed_offending_wells_.clear();
1153 const int episodeIdx = simulator_.episodeIndex();
1154 const auto& network = this->schedule()[episodeIdx].network();
1155 if (!this->wellsActive() && !network.active()) {
1160 if (iterationIdx == 0 && this->wellsActive()) {
1161 OPM_TIMEBLOCK(firstIterationAssmble);
1168 calculateExplicitQuantities(local_deferredLogger);
1169 prepareTimeStep(local_deferredLogger);
1172 "assemble() failed (It=0): ",
1173 this->terminal_output_, grid().comm());
1176 const bool well_group_control_changed = updateWellControlsAndNetwork(
false, dt, local_deferredLogger);
1180 if ( ! this->wellsActive() ) {
1184 assembleWellEqWithoutIteration(dt, local_deferredLogger);
1188 last_report_.well_group_control_changed = well_group_control_changed;
1189 last_report_.assemble_time_well += perfTimer.stop();
1195 template<
typename TypeTag>
1204 bool do_network_update =
true;
1205 bool well_group_control_changed =
false;
1206 Scalar network_imbalance = 0.0;
1208 const std::size_t iteration_to_relax = param_.network_max_strict_outer_iterations_;
1210 const std::size_t max_iteration = param_.network_max_outer_iterations_;
1211 std::size_t network_update_iteration = 0;
1212 while (do_network_update) {
1213 if (network_update_iteration >= max_iteration ) {
1215 const int episodeIdx = simulator_.episodeIndex();
1216 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1217 if (this->shouldBalanceNetwork(episodeIdx, iterationIdx + 1)) {
1218 const std::string msg = fmt::format(
"Maximum of {:d} network iterations has been used and we stop the update, \n"
1219 "and try again after the next Newton iteration (imbalance = {:.2e} bar, ctrl_change = {})",
1220 max_iteration, network_imbalance*1.0e-5, well_group_control_changed);
1221 local_deferredLogger.
debug(msg);
1223 if (this->terminal_output_) {
1224 const std::string msg = fmt::format(
"Maximum of {:d} network iterations has been used and we stop the update. \n"
1225 "The simulator will continue with unconverged network results (imbalance = {:.2e} bar, ctrl_change = {})",
1226 max_iteration, network_imbalance*1.0e-5, well_group_control_changed);
1227 local_deferredLogger.
info(msg);
1232 if (this->terminal_output_ && (network_update_iteration == iteration_to_relax) ) {
1233 local_deferredLogger.
debug(
"We begin using relaxed tolerance for network update now after " +
std::to_string(iteration_to_relax) +
" iterations ");
1235 const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
1237 const bool optimize_gas_lift = ( (network_update_iteration + 1) < std::max(max_iteration,
static_cast<std::size_t
>(2)) );
1238 std::tie(well_group_control_changed, do_network_update, network_imbalance) =
1239 updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, optimize_gas_lift, dt,local_deferredLogger);
1240 ++network_update_iteration;
1242 return well_group_control_changed;
1248 template<
typename TypeTag>
1249 std::tuple<bool, bool, typename BlackoilWellModel<TypeTag>::Scalar>
1252 const bool relax_network_tolerance,
1253 const bool optimize_gas_lift,
1258 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1259 const int reportStepIdx = simulator_.episodeIndex();
1260 this->updateAndCommunicateGroupData(reportStepIdx, iterationIdx,
1261 param_.nupcol_group_rate_tolerance_,
true, local_deferredLogger);
1262 const auto [more_inner_network_update, network_imbalance] =
1263 updateNetworks(mandatory_network_balance,
1264 local_deferredLogger,
1265 relax_network_tolerance);
1267 bool well_group_control_changed = updateWellControls(local_deferredLogger);
1269 bool alq_updated =
false;
1272 if (optimize_gas_lift) alq_updated = gaslift_.maybeDoGasLiftOptimize(simulator_,
1276 local_deferredLogger);
1278 prepareWellsBeforeAssembling(dt, local_deferredLogger);
1281 "updateWellControlsAndNetworkIteration() failed: ",
1282 this->terminal_output_, grid().comm());
1286 guideRateUpdateIsNeeded(reportStepIdx)) {
1287 const double simulationTime = simulator_.time();
1291 this->guide_rate_handler_.updateGuideRates(
1292 reportStepIdx, simulationTime, this->wellState(), this->groupState()
1297 const bool more_network_update = this->shouldBalanceNetwork(reportStepIdx, iterationIdx) &&
1298 (more_inner_network_update || well_group_control_changed || alq_updated);
1299 return {well_group_control_changed, more_network_update, network_imbalance};
1305 template <
typename TypeTag>
1311 const int reportStepIdx = this->simulator_.episodeIndex();
1312 const auto& network = this->schedule()[reportStepIdx].network();
1313 const auto& balance = this->schedule()[reportStepIdx].network_balance();
1314 const Scalar thp_tolerance = balance.thp_tolerance();
1316 if (!network.active()) {
1320 auto& well_state = this->wellState();
1321 auto& group_state = this->groupState();
1323 bool well_group_thp_updated =
false;
1324 for (
const std::string& nodeName : network.node_names()) {
1325 const bool has_choke = network.node(nodeName).as_choke();
1327 const auto& summary_state = this->simulator_.vanguard().summaryState();
1328 const Group& group = this->schedule().getGroup(nodeName, reportStepIdx);
1330 const auto pu = this->phase_usage_;
1332 std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
1333 Scalar gratTargetFromSales = 0.0;
1334 if (group_state.has_grat_sales_target(group.name()))
1335 gratTargetFromSales = group_state.grat_sales_target(group.name());
1337 const auto ctrl = group.productionControls(summary_state);
1338 auto cmode_tmp = ctrl.cmode;
1340 bool fld_none =
false;
1345 const Scalar efficiencyFactor = 1.0;
1346 const Group& parentGroup = this->schedule().getGroup(group.parent(), reportStepIdx);
1359 local_deferredLogger);
1360 target_tmp = target.first;
1361 cmode_tmp = target.second;
1363 const auto cmode = cmode_tmp;
1365 gratTargetFromSales, nodeName, group_state,
1366 group.has_gpmaint_control(cmode));
1370 target_tmp = tcalc.
groupTarget(ctrl, local_deferredLogger);
1373 const Scalar orig_target = target_tmp;
1375 auto mismatch = [&] (
auto group_thp) {
1378 for (
auto& well : this->well_container_) {
1379 std::string well_name = well->name();
1380 auto& ws = well_state.well(well_name);
1381 if (group.hasWell(well_name)) {
1382 well->setDynamicThpLimit(group_thp);
1383 const Well& well_ecl = this->wells_ecl_[well->indexOfWell()];
1384 const auto inj_controls = Well::InjectionControls(0);
1385 const auto prod_controls = well_ecl.productionControls(summary_state);
1386 well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger,
false,
false);
1391 return (group_rate - orig_target)/orig_target;
1394 const auto upbranch = network.uptree_branch(nodeName);
1395 const auto it = this->node_pressures_.find((*upbranch).uptree_node());
1396 const Scalar nodal_pressure = it->second;
1397 Scalar well_group_thp = nodal_pressure;
1399 std::optional<Scalar> autochoke_thp;
1400 if (
auto iter = this->well_group_thp_calc_.find(nodeName); iter != this->well_group_thp_calc_.end()) {
1401 autochoke_thp = this->well_group_thp_calc_.at(nodeName);
1405 std::array<Scalar, 2> range_initial;
1406 if (!autochoke_thp.has_value()){
1409 std::string node_name = nodeName;
1410 while (!network.node(node_name).terminal_pressure().has_value()) {
1411 auto branch = network.uptree_branch(node_name).value();
1412 node_name = branch.uptree_node();
1414 min_thp = network.node(node_name).terminal_pressure().value();
1418 std::array<Scalar, 2> range = {
Scalar{0.9}*min_thp,
Scalar{1.1}*max_thp};
1419 std::optional<Scalar> appr_sol;
1423 range_initial = {min_thp, max_thp};
1426 if (!autochoke_thp.has_value() || autochoke_thp.value() > nodal_pressure) {
1428 std::array<Scalar, 2> range = autochoke_thp.has_value() ?
1429 std::array<Scalar, 2>{
Scalar{0.9} * autochoke_thp.value(),
1430 Scalar{1.1} * autochoke_thp.value()} : range_initial;
1432 std::optional<Scalar> approximate_solution;
1433 const Scalar tolerance1 = thp_tolerance;
1434 local_deferredLogger.
debug(
"Using brute force search to bracket the group THP");
1437 if (approximate_solution.has_value()) {
1438 autochoke_thp = *approximate_solution;
1439 local_deferredLogger.
debug(
"Approximate group THP value found: " +
std::to_string(autochoke_thp.value()));
1440 }
else if (finding_bracket) {
1441 const Scalar tolerance2 = thp_tolerance;
1442 const int max_iteration_solve = 100;
1444 autochoke_thp = RegulaFalsiBisection<ThrowOnError>::
1445 solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration);
1450 autochoke_thp.reset();
1451 local_deferredLogger.
debug(
"Group THP solve failed due to bracketing failure");
1454 if (autochoke_thp.has_value()) {
1455 well_group_thp_calc_[nodeName] = autochoke_thp.value();
1458 well_group_thp = std::max(autochoke_thp.value(), nodal_pressure);
1461 for (
auto& well : this->well_container_) {
1462 std::string well_name = well->name();
1464 if (well->isInjector() || !well->wellEcl().predictionMode())
1467 if (group.hasWell(well_name)) {
1468 well->setDynamicThpLimit(well_group_thp);
1470 const auto& ws = this->wellState().well(well->indexOfWell());
1471 const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP;
1473 well->prepareWellBeforeAssembling(this->simulator_, dt, this->wellState(), this->groupState(), local_deferredLogger);
1478 const auto& current_well_group_thp = group_state.is_autochoke_group(nodeName) ? group_state.well_group_thp(nodeName) : 1e30;
1479 if (std::abs(current_well_group_thp - well_group_thp) > balance.pressure_tolerance()) {
1480 well_group_thp_updated =
true;
1481 group_state.update_well_group_thp(nodeName, well_group_thp);
1485 return well_group_thp_updated;
1488 template<
typename TypeTag>
1494 for (
auto& well : well_container_) {
1495 well->assembleWellEq(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1500 template<
typename TypeTag>
1506 for (
auto& well : well_container_) {
1507 well->prepareWellBeforeAssembling(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1512 template<
typename TypeTag>
1522 for (
auto& well: well_container_) {
1523 well->assembleWellEqWithoutIteration(simulator_, dt, this->wellState(), this->groupState(),
1527 this->terminal_output_, grid().comm());
1531#if COMPILE_GPU_BRIDGE
1532 template<
typename TypeTag>
1540 for(
unsigned int i = 0; i < well_container_.size(); i++){
1541 auto& well = well_container_[i];
1542 std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1544 wellContribs.
addNumBlocks(derived->linSys().getNumBlocks());
1549 wellContribs.
alloc();
1551 for(
unsigned int i = 0; i < well_container_.size(); i++){
1552 auto& well = well_container_[i];
1554 auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag>>(well);
1556 derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs);
1558 auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
1560 derived_ms->linSys().extract(wellContribs);
1562 OpmLog::warning(
"Warning unknown type of well");
1569 template<
typename TypeTag>
1574 for (
const auto& well: well_container_ ) {
1579 template<
typename TypeTag>
1584 const bool use_well_weights)
const
1586 int nw = this->numLocalWellsEnd();
1587 int rdofs = local_num_cells_;
1588 for (
int i = 0; i < nw; i++ ) {
1589 int wdof = rdofs + i;
1590 jacobian[wdof][wdof] = 1.0;
1593 for (
const auto& well : well_container_) {
1594 well->addWellPressureEquations(jacobian,
1602 template <
typename TypeTag>
1605 const std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress)
const
1610 for (
const auto& well : well_container_) {
1611 if (!well->isOperableAndSolvable() && !well->wellIsStopped()) {
1614 const auto& cells = well->cells();
1615 const auto& rates = well->connectionRates();
1616 for (
unsigned perfIdx = 0; perfIdx < rates.size(); ++perfIdx) {
1617 unsigned cellIdx = cells[perfIdx];
1618 auto rate = rates[perfIdx];
1621 using MatrixBlockType =
typename SparseMatrixAdapter::MatrixBlock;
1622 MatrixBlockType bMat(0.0);
1623 simulator_.model().linearizer().setResAndJacobi(res, bMat, rate);
1624 residual[cellIdx] += res;
1625 *diagMatAddress[cellIdx] += bMat;
1631 template<
typename TypeTag>
1636 int nw = this->numLocalWellsEnd();
1637 int rdofs = local_num_cells_;
1638 for (
int i = 0; i < nw; ++i) {
1639 int wdof = rdofs + i;
1640 jacobian.entry(wdof,wdof) = 1.0;
1642 const auto wellconnections = this->getMaxWellConnections();
1643 for (
int i = 0; i < nw; ++i) {
1644 const auto& perfcells = wellconnections[i];
1645 for (
int perfcell : perfcells) {
1646 int wdof = rdofs + i;
1647 jacobian.entry(wdof, perfcell) = 0.0;
1648 jacobian.entry(perfcell, wdof) = 0.0;
1654 template<
typename TypeTag>
1662 for (
const auto& well : well_container_) {
1663 const auto& cells = well->cells();
1664 x_local_.resize(cells.size());
1666 for (
size_t i = 0; i < cells.size(); ++i) {
1667 x_local_[i] = x[cells[i]];
1669 well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_,
1670 this->wellState(), local_deferredLogger);
1674 "recoverWellSolutionAndUpdateWellState() failed: ",
1675 this->terminal_output_, simulator_.vanguard().grid().comm());
1679 template<
typename TypeTag>
1685 OPM_THROW(std::logic_error,
"Attempt to call NLDD method without a NLDD solver");
1688 return nldd_->recoverWellSolutionAndUpdateWellState(x, domainIdx);
1692 template<
typename TypeTag>
1695 getWellConvergence(
const std::vector<Scalar>& B_avg,
bool checkWellGroupControls)
const
1701 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1702 for (
const auto& well : well_container_) {
1703 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
1704 local_report += well->getWellConvergence(
1705 simulator_, this->wellState(), B_avg, local_deferredLogger,
1706 iterationIdx > param_.strict_outer_iter_wells_);
1710 report.
setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
1711 local_report += report;
1720 if (checkWellGroupControls) {
1724 if (this->terminal_output_) {
1729 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
1730 OpmLog::debug(
"NaN residual found with phase " +
std::to_string(f.phase()) +
" for well " + f.wellName());
1731 }
else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
1732 OpmLog::debug(
"Too large residual found with phase " +
std::to_string(f.phase()) +
" for well " + f.wellName());
1743 template<
typename TypeTag>
1749 for (
auto& well : well_container_) {
1758 template<
typename TypeTag>
1764 if (!this->wellsActive()) {
1767 const int episodeIdx = simulator_.episodeIndex();
1768 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1769 const auto& comm = simulator_.vanguard().grid().comm();
1771 bool changed_well_group =
false;
1772 const Group& fieldGroup = this->schedule().getGroup(
"FIELD", episodeIdx);
1775 const std::size_t max_iter = param_.well_group_constraints_max_iterations_;
1776 while(!changed_well_group && iter < max_iter) {
1777 changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
1780 bool changed_well_to_group =
false;
1782 OPM_TIMEBLOCK(UpdateWellControls);
1786 for (
const auto& well : well_container_) {
1788 const bool changed_well = well->
updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1790 changed_well_to_group = changed_well || changed_well_to_group;
1794 simulator_.gridView().comm());
1797 changed_well_to_group = comm.sum(
static_cast<int>(changed_well_to_group));
1798 if (changed_well_to_group) {
1799 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1800 changed_well_group =
true;
1804 bool changed_well_individual =
false;
1809 for (
const auto& well : well_container_) {
1811 const bool changed_well = well->
updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1813 changed_well_individual = changed_well || changed_well_individual;
1817 simulator_.gridView().comm());
1820 changed_well_individual = comm.sum(
static_cast<int>(changed_well_individual));
1821 if (changed_well_individual) {
1822 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1823 changed_well_group =
true;
1829 this->updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState());
1831 return changed_well_group;
1835 template<
typename TypeTag>
1836 std::tuple<bool, typename BlackoilWellModel<TypeTag>::Scalar>
1840 const bool relax_network_tolerance)
1843 const int episodeIdx = simulator_.episodeIndex();
1844 const auto& network = this->schedule()[episodeIdx].network();
1845 if (!this->wellsActive() && !network.active()) {
1846 return {
false, 0.0};
1849 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1850 const auto& comm = simulator_.vanguard().grid().comm();
1853 Scalar network_imbalance = 0.0;
1854 bool more_network_update =
false;
1855 if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
1856 OPM_TIMEBLOCK(BalanceNetwork);
1857 const double dt = this->simulator_.timeStepSize();
1859 const bool well_group_thp_updated = computeWellGroupThp(dt, deferred_logger);
1860 const int max_number_of_sub_iterations = param_.network_max_sub_iterations_;
1861 const Scalar network_pressure_update_damping_factor = param_.network_pressure_update_damping_factor_;
1862 const Scalar network_max_pressure_update = param_.network_max_pressure_update_in_bars_ * unit::barsa;
1863 bool more_network_sub_update =
false;
1864 for (
int i = 0; i < max_number_of_sub_iterations; i++) {
1865 const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, network_pressure_update_damping_factor, network_max_pressure_update);
1866 network_imbalance = comm.max(local_network_imbalance);
1867 const auto& balance = this->schedule()[episodeIdx].network_balance();
1868 constexpr Scalar relaxation_factor = 10.0;
1869 const Scalar tolerance = relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() : balance.pressure_tolerance();
1870 more_network_sub_update = this->networkActive() && network_imbalance > tolerance;
1871 if (!more_network_sub_update)
1874 for (
const auto& well : well_container_) {
1875 if (well->isInjector() || !well->wellEcl().predictionMode())
1878 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
1879 if (it != this->node_pressures_.end()) {
1880 const auto& ws = this->wellState().well(well->indexOfWell());
1881 const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP;
1883 well->prepareWellBeforeAssembling(this->simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1887 this->updateAndCommunicateGroupData(episodeIdx, iterationIdx, param_.nupcol_group_rate_tolerance_,
1888 true, deferred_logger);
1890 more_network_update = more_network_sub_update || well_group_thp_updated;
1892 return { more_network_update, network_imbalance };
1896 template<
typename TypeTag>
1900 const int iterationIdx,
1903 this->updateAndCommunicateGroupData(reportStepIdx,
1905 param_.nupcol_group_rate_tolerance_,
1913 for (
const auto& well : well_container_) {
1915 const auto& ws = this->wellState().well(well->indexOfWell());
1916 if (ws.production_cmode == Well::ProducerCMode::GRUP ||
1917 ws.injection_cmode == Well::InjectorCMode::GRUP)
1919 well->updateWellStateWithTarget(simulator_, this->groupState(),
1920 this->wellState(), deferred_logger);
1924 simulator_.gridView().comm())
1925 this->updateAndCommunicateGroupData(reportStepIdx,
1927 param_.nupcol_group_rate_tolerance_,
1932 template<
typename TypeTag>
1937 const int reportStepIdx,
1938 const int iterationIdx)
1941 bool changed =
false;
1943 const int nupcol = this->schedule()[reportStepIdx].nupcol();
1944 const int max_number_of_group_switches = iterationIdx < nupcol ? 9999 : param_.max_number_of_group_switches_;
1945 bool changed_hc = this->checkGroupHigherConstraints( group, deferred_logger, reportStepIdx, max_number_of_group_switches);
1948 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1951 bool changed_individual =
1953 updateGroupIndividualControl(group,
1955 max_number_of_group_switches,
1956 this->switched_inj_groups_,
1957 this->switched_prod_groups_,
1958 this->closed_offending_wells_,
1963 if (changed_individual) {
1965 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1968 for (
const std::string& groupName : group.groups()) {
1969 bool changed_this = updateGroupControls( this->schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx,iterationIdx);
1970 changed = changed || changed_this;
1975 template<
typename TypeTag>
1982 for (
const auto& well : well_container_) {
1983 const auto& wname = well->name();
1984 const auto wasClosed = wellTestState.well_is_closed(wname);
1985 well->checkWellOperability(simulator_,
1987 local_deferredLogger);
1988 const bool under_zero_target =
1989 well->wellUnderZeroGroupRateTarget(this->simulator_,
1991 local_deferredLogger);
1992 well->updateWellTestState(this->wellState().well(wname),
1997 local_deferredLogger);
1999 if (!wasClosed && wellTestState.well_is_closed(wname)) {
2000 this->closed_this_step_.insert(wname);
2004 for (
const auto& [group_name, to] : this->closed_offending_wells_) {
2005 if (this->hasOpenLocalWell(to.second) &&
2006 !this->wasDynamicallyShutThisTimeStep(to.second))
2008 wellTestState.close_well(to.second,
2009 WellTestConfig::Reason::GROUP,
2011 this->updateClosedWellsThisStep(to.second);
2012 const std::string msg =
2013 fmt::format(
"Procedure on exceeding {} limit is WELL for group {}. "
2019 local_deferredLogger.
info(msg);
2027 if (this->terminal_output_) {
2033 template<
typename TypeTag>
2037 std::string& exc_msg,
2042 const int np = this->numPhases();
2043 std::vector<Scalar> potentials;
2044 const auto& well = well_container_[widx];
2045 std::string cur_exc_msg;
2048 well->computeWellPotentials(simulator_, well_state_copy, potentials, deferred_logger);
2053 exc_msg += fmt::format(
"\nFor well {}: {}", well->name(), cur_exc_msg);
2055 exc_type = std::max(exc_type, cur_exc_type);
2059 auto& ws = this->wellState().well(well->indexOfWell());
2060 for (
int p = 0; p < np; ++p) {
2062 ws.well_potentials[p] = std::max(
Scalar{0.0}, potentials[p]);
2068 template <
typename TypeTag>
2073 for (
const auto& wellPtr : this->well_container_) {
2074 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2082 template <
typename TypeTag>
2093 for (
const auto& shutWell : this->local_shut_wells_) {
2094 if (!this->wells_ecl_[shutWell].hasConnections()) {
2099 auto wellPtr = this->
template createTypedWellPointer
2102 wellPtr->
init(&this->phase_usage_, this->depth_, this->gravity_, this->B_avg_,
true);
2104 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2112 template <
typename TypeTag>
2126 template<
typename TypeTag>
2132 const auto episodeIdx = simulator_.episodeIndex();
2133 this->updateNetworkActiveState(episodeIdx);
2137 const bool do_prestep_network_rebalance = param_.pre_solve_network_ && this->needPreStepNetworkRebalance(episodeIdx);
2139 for (
const auto& well : well_container_) {
2140 auto& events = this->wellState().well(well->indexOfWell()).events;
2142 well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), deferred_logger);
2143 well->updatePrimaryVariables(simulator_, this->wellState(), deferred_logger);
2149 if (events.hasEvent(ScheduleEvents::REQUEST_OPEN_WELL)) {
2150 events.clearEvent(ScheduleEvents::REQUEST_OPEN_WELL);
2153 if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
2155 well->solveWellEquation(simulator_, this->wellState(), this->groupState(), deferred_logger);
2156 }
catch (
const std::exception& e) {
2157 const std::string msg =
"Compute initial well solution for " + well->name() +
" initially failed. Continue with the previous rates";
2158 deferred_logger.
warning(
"WELL_INITIAL_SOLVE_FAILED", msg);
2163 well->resetWellOperability();
2165 updatePrimaryVariables(deferred_logger);
2168 if (do_prestep_network_rebalance) doPreStepNetworkRebalance(deferred_logger);
2171 template<
typename TypeTag>
2176 std::vector< Scalar > B_avg(numComponents(),
Scalar() );
2177 const auto& grid = simulator_.vanguard().grid();
2178 const auto& gridView = grid.leafGridView();
2182 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2183 elemCtx.updatePrimaryStencil(elem);
2184 elemCtx.updatePrimaryIntensiveQuantities(0);
2186 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
2187 const auto& fs = intQuants.fluidState();
2189 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
2191 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2195 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2196 auto& B = B_avg[ compIdx ];
2198 B += 1 / fs.invB(phaseIdx).value();
2200 if constexpr (has_solvent_) {
2201 auto& B = B_avg[solventSaturationIdx];
2202 B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
2208 grid.comm().sum(B_avg.data(), B_avg.size());
2209 B_avg_.resize(B_avg.size());
2210 std::transform(B_avg.begin(), B_avg.end(), B_avg_.begin(),
2211 [gcells = global_num_cells_](
const auto bval)
2212 { return bval / gcells; });
2219 template<
typename TypeTag>
2224 for (
const auto& well : well_container_) {
2229 template<
typename TypeTag>
2233 const auto& grid = simulator_.vanguard().
grid();
2234 const auto& eclProblem = simulator_.problem();
2235 const unsigned numCells = grid.size(0);
2237 this->pvt_region_idx_.resize(numCells);
2238 for (
unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
2239 this->pvt_region_idx_[cellIdx] =
2240 eclProblem.pvtRegionIndex(cellIdx);
2245 template<
typename TypeTag>
2255 int numComp = this->numPhases() < 3 ? this->numPhases() : FluidSystem::numComponents;
2256 if constexpr (has_solvent_) {
2262 template<
typename TypeTag>
2266 const auto& eclProblem = simulator_.problem();
2267 depth_.resize(local_num_cells_);
2268 for (
unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
2269 depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
2273 template<
typename TypeTag>
2276 getWell(
const std::string& well_name)
const
2279 auto well = std::find_if(well_container_.begin(),
2280 well_container_.end(),
2282 return elem->name() == well_name;
2285 assert(well != well_container_.end());
2290 template <
typename TypeTag>
2295 return std::max(this->simulator_.episodeIndex(), 0);
2302 template<
typename TypeTag>
2307 const std::vector<Scalar>& production_rates,
2308 std::vector<Scalar>& resv_coeff)
2310 rateConverter_->calcCoeff(fipnum, pvtreg, production_rates, resv_coeff);
2313 template<
typename TypeTag>
2318 std::vector<Scalar>& resv_coeff)
2320 rateConverter_->calcInjCoeff(fipnum, pvtreg, resv_coeff);
2324 template <
typename TypeTag>
2329 if constexpr (has_energy_) {
2330 int np = this->numPhases();
2331 Scalar cellInternalEnergy;
2335 const int nw = this->numLocalWells();
2336 for (
auto wellID = 0*nw; wellID < nw; ++wellID) {
2337 const Well& well = this->wells_ecl_[wellID];
2338 auto& ws = this->wellState().well(wellID);
2339 if (well.isInjector()) {
2340 if (ws.status != WellStatus::STOP) {
2341 this->wellState().well(wellID).temperature = well.inj_temperature();
2346 std::array<Scalar,2> weighted{0.0,0.0};
2347 auto& [weighted_temperature, total_weight] = weighted;
2349 auto& well_info = this->local_parallel_well_info_[wellID].get();
2350 auto& perf_data = ws.perf_data;
2351 auto& perf_phase_rate = perf_data.phase_rates;
2353 using int_type =
decltype(this->well_perf_data_[wellID].size());
2354 for (int_type perf = 0, end_perf = this->well_perf_data_[wellID].size(); perf < end_perf; ++perf) {
2355 const int cell_idx = this->well_perf_data_[wellID][perf].cell_index;
2356 const auto& intQuants = simulator_.model().intensiveQuantities(cell_idx, 0);
2357 const auto& fs = intQuants.fluidState();
2360 Scalar cellTemperatures = fs.temperature(0).value();
2362 Scalar weight_factor = 0.0;
2363 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2364 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2367 cellInternalEnergy = fs.enthalpy(phaseIdx).value() -
2368 fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
2369 cellBinv = fs.invB(phaseIdx).value();
2370 cellDensity = fs.density(phaseIdx).value();
2371 perfPhaseRate = perf_phase_rate[perf*np + phaseIdx];
2372 weight_factor += cellDensity * perfPhaseRate / cellBinv * cellInternalEnergy / cellTemperatures;
2374 weight_factor = std::abs(weight_factor) + 1e-13;
2375 total_weight += weight_factor;
2376 weighted_temperature += weight_factor * cellTemperatures;
2378 well_info.communication().sum(weighted.data(), 2);
2379 this->wellState().well(wellID).temperature = weighted_temperature / total_weight;
2385 template <
typename TypeTag>
2389 const auto reportStepIdx =
static_cast<unsigned int>(this->reportStepIndex());
2390 const auto& trMod = this->simulator_.problem().tracerModel();
2396 this->assignMswTracerRates(wsrpt, trMod.getMswTracerRates(), reportStepIdx);
#define OPM_END_PARALLEL_TRY_CATCH_LOG(obptc_logger, obptc_prefix, obptc_output, comm)
Catch exception, log, and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:202
#define OPM_DEFLOG_THROW(Exception, message, deferred_logger)
Definition: DeferredLoggingErrorHelpers.hpp:45
#define OPM_END_PARALLEL_TRY_CATCH(prefix, comm)
Catch exception and throw in a parallel try-catch clause.
Definition: DeferredLoggingErrorHelpers.hpp:192
#define OPM_PARALLEL_CATCH_CLAUSE(obptc_exc_type, obptc_exc_msg)
Inserts catch classes for the parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:166
#define OPM_BEGIN_PARALLEL_TRY_CATCH()
Macro to setup the try of a parallel try-catch.
Definition: DeferredLoggingErrorHelpers.hpp:158
void logAndCheckForExceptionsAndThrow(Opm::DeferredLogger &deferred_logger, Opm::ExceptionType::ExcEnum exc_type, const std::string &message, const bool terminal_output, Opm::Parallel::Communication comm)
Definition: DeferredLoggingErrorHelpers.hpp:111
Class for handling constraints for the blackoil well model.
Definition: BlackoilWellModelConstraints.hpp:41
Class for handling the gaslift in the blackoil well model.
Definition: BlackoilWellModelGasLift.hpp:94
Class for handling the blackoil well model.
Definition: BlackoilWellModelGeneric.hpp:94
Class for handling the guide rates in the blackoil well model.
Definition: BlackoilWellModelGuideRates.hpp:46
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:102
void initializeGroupStructure(const int reportStepIdx)
Definition: BlackoilWellModel_impl.hpp:298
void init()
Definition: BlackoilWellModel_impl.hpp:165
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: BlackoilWellModel.hpp:111
void initWellContainer(const int reportStepIdx) override
Definition: BlackoilWellModel_impl.hpp:184
void beginReportStep(const int time_step)
Definition: BlackoilWellModel_impl.hpp:201
Dune::FieldVector< Scalar, numEq > VectorBlockType
Definition: BlackoilWellModel.hpp:130
GetPropType< TypeTag, Properties::ElementContext > ElementContext
Definition: BlackoilWellModel.hpp:108
GetPropType< TypeTag, Properties::Grid > Grid
Definition: BlackoilWellModel.hpp:105
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: BlackoilWellModel.hpp:110
void initializeWellState(const int timeStepIdx)
Definition: BlackoilWellModel_impl.hpp:761
const Grid & grid() const
Definition: BlackoilWellModel.hpp:367
const SimulatorReportSingle & lastReport() const
Definition: BlackoilWellModel_impl.hpp:635
void addWellContributions(SparseMatrixAdapter &jacobian) const
Definition: BlackoilWellModel_impl.hpp:1572
Dune::BCRSMatrix< Opm::MatrixBlock< Scalar, 1, 1 > > PressureMatrix
Definition: BlackoilWellModel.hpp:295
bool empty() const
Definition: BlackoilWellModel.hpp:340
void computeTotalRatesForDof(RateVector &rate, unsigned globalIdx) const
Definition: BlackoilWellModel_impl.hpp:722
void beginTimeStep()
Definition: BlackoilWellModel_impl.hpp:340
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: BlackoilWellModel.hpp:112
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Definition: BlackoilWellModel_impl.hpp:1746
void updatePrimaryVariables(DeferredLogger &deferred_logger)
Definition: BlackoilWellModel_impl.hpp:2222
void initializeLocalWellStructure(const int reportStepIdx, const bool enableWellPIScaling)
Definition: BlackoilWellModel_impl.hpp:253
Dune::BlockVector< VectorBlockType > BVector
Definition: BlackoilWellModel.hpp:131
BlackoilWellModel(Simulator &simulator)
Definition: BlackoilWellModel_impl.hpp:157
void wellTesting(const int timeStepIdx, const double simulationTime, DeferredLogger &deferred_logger)
Definition: BlackoilWellModel_impl.hpp:545
void timeStepSucceeded(const double simulationTime, const double dt)
Definition: BlackoilWellModel_impl.hpp:645
std::shared_ptr< WellInterface< TypeTag > > WellInterfacePtr
Definition: BlackoilWellModel.hpp:188
void createWellContainer(const int report_step) override
Definition: BlackoilWellModel_impl.hpp:804
void endReportStep()
Definition: BlackoilWellModel_impl.hpp:618
Definition: ConvergenceReport.hpp:38
void setWellFailed(const WellFailure &wf)
Definition: ConvergenceReport.hpp:270
void setWellGroupTargetsViolated(const bool wellGroupTargetsViolated)
Definition: ConvergenceReport.hpp:288
const std::vector< WellFailure > & wellFailures() const
Definition: ConvergenceReport.hpp:369
Definition: DeferredLogger.hpp:57
void info(const std::string &tag, const std::string &message)
void warning(const std::string &tag, const std::string &message)
void debug(const std::string &tag, const std::string &message)
std::map< std::string, std::pair< const Well *, int > > GLiftEclWells
Definition: GasLiftGroupInfo.hpp:70
Definition: StandardWell.hpp:60
virtual void init(const PhaseUsage *phase_usage_arg, const std::vector< Scalar > &depth_arg, const Scalar gravity_arg, const std::vector< Scalar > &B_avg, const bool changed_to_open_this_step) override
Definition: StandardWell_impl.hpp:76
Definition: TargetCalculator.hpp:43
RateType calcModeRateFromRates(const std::vector< RateType > &rates) const
Definition: TargetCalculator.hpp:54
Scalar groupTarget(const std::optional< Group::ProductionControls > &ctrl, DeferredLogger &deferred_logger) const
Class for computing BHP limits.
Definition: WellBhpThpCalculator.hpp:41
Definition: WellContributions.hpp:51
void alloc()
Allocate memory for the StandardWells.
void setBlockSize(unsigned int dim, unsigned int dim_wells)
void addNumBlocks(unsigned int numBlocks)
Class for computing well group controls.
Definition: WellGroupControls.hpp:49
static void setCmodeGroup(const Group &group, const Schedule &schedule, const SummaryState &summaryState, const int reportStepIdx, GroupState< Scalar > &group_state)
static void setRegionAveragePressureCalculator(const Group &group, const Schedule &schedule, const int reportStepIdx, const FieldPropsManager &fp, const PhaseUsage &pu, std::map< std::string, std::unique_ptr< AverageRegionalPressureType > > ®ionalAveragePressureCalculator)
static void updateGpMaintTargetForGroups(const Group &group, const Schedule &schedule, const RegionalValues ®ional_values, const int reportStepIdx, const double dt, const WellState< Scalar > &well_state, GroupState< Scalar > &group_state)
static void accumulateGroupEfficiencyFactor(const Group &group, const Schedule &schedule, const int reportStepIdx, Scalar &factor)
int indexOfWell() const
Index of well in the wells struct and wellState.
Definition: WellInterface.hpp:77
bool updateWellControl(const Simulator &simulator, const IndividualOrGroup iog, WellState< Scalar > &well_state, const GroupState< Scalar > &group_state, DeferredLogger &deferred_logger)
Definition: WellInterface_impl.hpp:191
virtual void updateProductivityIndex(const Simulator &simulator, const WellProdIndexCalculator< Scalar > &wellPICalc, WellState< Scalar > &well_state, DeferredLogger &deferred_logger) const =0
Definition: WellState.hpp:65
ExcEnum
Definition: DeferredLogger.hpp:45
@ NONE
Definition: DeferredLogger.hpp:46
Dune::Communication< MPIComm > Communication
Definition: ParallelCommunication.hpp:30
Definition: blackoilboundaryratevector.hh:39
Opm::DeferredLogger gatherDeferredLogger(const Opm::DeferredLogger &local_deferredlogger, Parallel::Communication communicator)
Create a global log combining local logs.
ConvergenceReport gatherConvergenceReport(const ConvergenceReport &local_report, Parallel::Communication communicator)
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
A struct for returning timing data from a simulator to its caller.
Definition: SimulatorReport.hpp:34