24#ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
25#define OPM_BLACKOILWELLMODEL_IMPL_HEADER_INCLUDED
27#include <opm/simulators/wells/BlackoilWellModel.hpp>
30#include <opm/grid/utility/cartesianToCompressed.hpp>
31#include <opm/common/utility/numeric/RootFinders.hpp>
33#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
34#include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
35#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
36#include <opm/input/eclipse/Schedule/Well/WellTestConfig.hpp>
38#include <opm/input/eclipse/Units/UnitSystem.hpp>
40#include <opm/simulators/wells/BlackoilWellModelConstraints.hpp>
41#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
42#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
43#include <opm/simulators/wells/VFPProperties.hpp>
44#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
45#include <opm/simulators/wells/WellGroupHelpers.hpp>
46#include <opm/simulators/wells/TargetCalculator.hpp>
48#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
49#include <opm/simulators/utils/MPIPacker.hpp>
50#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
53#include <opm/simulators/linalg/bda/WellContributions.hpp>
57#include <opm/simulators/utils/MPISerializer.hpp>
66#include <fmt/format.h>
69 template<
typename TypeTag>
70 BlackoilWellModel<TypeTag>::
71 BlackoilWellModel(Simulator& simulator,
const PhaseUsage& phase_usage)
72 : BlackoilWellModelGeneric<Scalar>(simulator.vanguard().schedule(),
73 simulator.vanguard().summaryState(),
74 simulator.vanguard().eclState(),
76 simulator.gridView().comm())
77 , simulator_(simulator)
79 this->terminal_output_ = (simulator.gridView().comm().rank() == 0)
80 && Parameters::Get<Parameters::EnableTerminalOutput>();
82 local_num_cells_ = simulator_.gridView().size(0);
85 global_num_cells_ = simulator_.vanguard().globalNumCells();
88 auto& parallel_wells = simulator.vanguard().parallelWells();
90 this->parallel_well_info_.reserve(parallel_wells.size());
91 for(
const auto& name_bool : parallel_wells) {
92 this->parallel_well_info_.emplace_back(name_bool, grid().comm());
96 this->alternative_well_rate_init_ =
97 Parameters::Get<Parameters::AlternativeWellRateInit>();
99 using SourceDataSpan =
100 typename PAvgDynamicSourceData<Scalar>::template SourceDataSpan<Scalar>;
102 this->wbpCalculationService_
103 .localCellIndex([
this](
const std::size_t globalIndex)
104 {
return this->compressedIndexForInterior(globalIndex); })
105 .evalCellSource([
this](
const int localCell, SourceDataSpan sourceTerms)
107 using Item =
typename SourceDataSpan::Item;
109 const auto* intQuants = this->simulator_.model()
110 .cachedIntensiveQuantities(localCell, 0);
111 const auto& fs = intQuants->fluidState();
113 sourceTerms.set(Item::PoreVol, intQuants->porosity().value() *
114 this->simulator_.model().dofTotalVolume(localCell));
116 constexpr auto io = FluidSystem::oilPhaseIdx;
117 constexpr auto ig = FluidSystem::gasPhaseIdx;
118 constexpr auto iw = FluidSystem::waterPhaseIdx;
121 const auto haveOil = FluidSystem::phaseIsActive(io);
122 const auto haveGas = FluidSystem::phaseIsActive(ig);
123 const auto haveWat = FluidSystem::phaseIsActive(iw);
125 auto weightedPhaseDensity = [&fs](
const auto ip)
127 return fs.saturation(ip).value() * fs.density(ip).value();
130 if (haveOil) { sourceTerms.set(Item::Pressure, fs.pressure(io).value()); }
131 else if (haveGas) { sourceTerms.set(Item::Pressure, fs.pressure(ig).value()); }
132 else { sourceTerms.set(Item::Pressure, fs.pressure(iw).value()); }
136 if (haveOil) { rho += weightedPhaseDensity(io); }
137 if (haveGas) { rho += weightedPhaseDensity(ig); }
138 if (haveWat) { rho += weightedPhaseDensity(iw); }
140 sourceTerms.set(Item::MixtureDensity, rho);
144 template<
typename TypeTag>
145 BlackoilWellModel<TypeTag>::
146 BlackoilWellModel(Simulator& simulator) :
151 template<
typename TypeTag>
153 BlackoilWellModel<TypeTag>::
156 extractLegacyCellPvtRegionIndex_();
157 extractLegacyDepth_();
159 gravity_ = simulator_.problem().gravity()[2];
161 this->initial_step_ =
true;
164 simulator_.model().addAuxiliaryModule(
this);
166 is_cell_perforated_.resize(local_num_cells_,
false);
170 template<
typename TypeTag>
172 BlackoilWellModel<TypeTag>::
173 initWellContainer(
const int reportStepIdx)
175 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
176 + ScheduleEvents::NEW_WELL;
177 const auto& events = this->schedule()[reportStepIdx].wellgroup_events();
178 for (
auto& wellPtr : this->well_container_) {
179 const bool well_opened_this_step = this->report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask);
180 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_, this->B_avg_, well_opened_this_step);
185 template<
typename TypeTag>
190 if (!param_.matrix_add_well_contributions_) {
195 const auto& schedule_wells = this->schedule().getWellsatEnd();
196 auto possibleFutureConnections = this->schedule().getPossibleFutureConnections();
200 const auto& comm = this->simulator_.vanguard().grid().comm();
202 ser.
broadcast(possibleFutureConnections);
205 for (
const auto& well : schedule_wells)
207 std::vector<int> wellCells = this->getCellsForConnections(well);
209 const auto possibleFutureConnectionSetIt = possibleFutureConnections.find(well.name());
210 if (possibleFutureConnectionSetIt != possibleFutureConnections.end()) {
211 for (
auto& global_index : possibleFutureConnectionSetIt->second) {
212 int compressed_idx = compressedIndexForInterior(global_index);
213 if (compressed_idx >= 0) {
214 wellCells.push_back(compressed_idx);
218 for (
int cellIdx : wellCells) {
219 neighbors[cellIdx].insert(wellCells.begin(),
226 template<
typename TypeTag>
229 linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
231 OPM_BEGIN_PARALLEL_TRY_CATCH();
232 for (
const auto& well: well_container_) {
235 if (param_.matrix_add_well_contributions_) {
236 well->addWellContributions(jacobian);
241 const auto& cells = well->cells();
242 linearize_res_local_.resize(cells.size());
244 for (
size_t i = 0; i < cells.size(); ++i) {
245 linearize_res_local_[i] = res[cells[i]];
248 well->apply(linearize_res_local_);
250 for (
size_t i = 0; i < cells.size(); ++i) {
251 res[cells[i]] = linearize_res_local_[i];
254 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel::linearize failed: ",
255 simulator_.gridView().comm());
259 template<
typename TypeTag>
262 linearizeDomain(
const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res)
267 for (
const auto& well: well_container_) {
268 if (well_domain_.at(well->name()) == domain.index) {
271 if (param_.matrix_add_well_contributions_) {
272 well->addWellContributions(jacobian);
277 const auto& cells = well->cells();
278 linearize_res_local_.resize(cells.size());
280 for (
size_t i = 0; i < cells.size(); ++i) {
281 linearize_res_local_[i] = res[cells[i]];
284 well->apply(linearize_res_local_);
286 for (
size_t i = 0; i < cells.size(); ++i) {
287 res[cells[i]] = linearize_res_local_[i];
294 template<
typename TypeTag>
296 BlackoilWellModel<TypeTag>::
297 beginReportStep(
const int timeStepIdx)
299 DeferredLogger local_deferredLogger{};
301 this->report_step_starts_ =
true;
304 this->rateConverter_ = std::make_unique<RateConverterType>
305 (this->phase_usage_, std::vector<int>(this->local_num_cells_, 0));
309 const auto enableWellPIScaling =
true;
310 this->initializeLocalWellStructure(timeStepIdx, enableWellPIScaling);
313 this->initializeGroupStructure(timeStepIdx);
315 const auto& comm = this->simulator_.vanguard().grid().comm();
317 OPM_BEGIN_PARALLEL_TRY_CATCH()
321 this->rateConverter_->template defineState<ElementContext>(this->simulator_);
325 const auto& sched_state = this->schedule()[timeStepIdx];
327 this->vfp_properties_ = std::make_unique<VFPProperties<Scalar>>
328 (sched_state.vfpinj(), sched_state.vfpprod(), this->wellState());
331 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
332 "beginReportStep() failed: ",
333 this->terminal_output_, comm)
337 this->commitWGState();
339 this->wellStructureChangedDynamically_ =
false;
346 template <
typename TypeTag>
350 const bool enableWellPIScaling)
354 const auto& comm = this->simulator_.vanguard().grid().comm();
357 this->wells_ecl_ = this->getLocalWells(reportStepIdx);
358 this->local_parallel_well_info_ =
359 this->createLocalParallelWellInfo(this->wells_ecl_);
364 OPM_BEGIN_PARALLEL_TRY_CATCH()
366 this->initializeWellPerfData();
367 this->initializeWellState(reportStepIdx);
368 this->initializeWBPCalculationService();
370 if (this->param_.use_multisegment_well_ && this->anyMSWellOpenLocal()) {
371 this->wellState().initWellStateMSWell(this->wells_ecl_, &this->prevWellState());
374 this->initializeWellProdIndCalculators();
376 if (enableWellPIScaling && this->schedule()[reportStepIdx].events()
377 .hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX))
379 this->runWellPIScaling(reportStepIdx, local_deferredLogger);
382 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
383 "Failed to initialize local well structure: ",
384 this->terminal_output_, comm)
391 template <
typename TypeTag>
398 const auto& comm = this->simulator_.vanguard().grid().comm();
400 OPM_BEGIN_PARALLEL_TRY_CATCH()
402 const auto& fieldGroup =
403 this->schedule().getGroup(
"FIELD", reportStepIdx);
407 this->summaryState(),
413 if (this->schedule()[reportStepIdx].has_gpmaint()) {
418 this->eclState_.fieldProps(),
420 this->regionalAveragePressureCalculator_);
423 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
424 "Failed to initialize group structure: ",
425 this->terminal_output_, comm)
433 template<
typename TypeTag>
438 OPM_TIMEBLOCK(beginTimeStep);
440 this->updateAverageFormationFactor();
444 this->switched_prod_groups_.clear();
445 this->switched_inj_groups_.clear();
447 if (this->wellStructureChangedDynamically_) {
452 const auto reportStepIdx =
453 this->simulator_.episodeIndex();
457 const auto enableWellPIScaling =
false;
459 this->initializeLocalWellStructure(reportStepIdx, enableWellPIScaling);
460 this->initializeGroupStructure(reportStepIdx);
462 this->commitWGState();
468 this->wellStructureChangedDynamically_ =
false;
471 this->resetWGState();
473 const int reportStepIdx = simulator_.episodeIndex();
474 this->updateAndCommunicateGroupData(reportStepIdx,
475 simulator_.model().newtonMethod().numIterations());
477 this->wellState().updateWellsDefaultALQ(this->schedule(), reportStepIdx, this->summaryState());
478 this->wellState().gliftTimeStepInit();
480 const double simulationTime = simulator_.time();
481 OPM_BEGIN_PARALLEL_TRY_CATCH();
484 wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
487 createWellContainer(reportStepIdx);
490 const Grid& grid = simulator_.vanguard().grid();
491 this->wells_active_ = grid.comm().max(!this->well_container_.empty());
496 this->initWellContainer(reportStepIdx);
499 std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(),
false);
500 for (
auto& well : well_container_) {
501 well->updatePerforatedCell(is_cell_perforated_);
505 this->calculateEfficiencyFactors(reportStepIdx);
507 if constexpr (has_polymer_)
509 if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
510 this->setRepRadiusPerfLength();
515 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
"beginTimeStep() failed: ",
516 this->terminal_output_, simulator_.vanguard().grid().comm());
518 for (
auto& well : well_container_) {
519 well->setVFPProperties(this->vfp_properties_.get());
520 well->setGuideRate(&this->guideRate_);
523 this->updateFiltrationModelsPreStep(local_deferredLogger);
526 for (
auto& well : well_container_) {
527 well->closeCompletions(this->wellTestState());
533 const auto& summaryState = simulator_.vanguard().summaryState();
534 if (alternative_well_rate_init_) {
539 for (
auto& well : well_container_) {
540 const bool zero_target = well->stoppedOrZeroRateTarget(simulator_, this->wellState(), local_deferredLogger);
541 if (well->isProducer() && !zero_target) {
542 well->updateWellStateRates(simulator_, this->wellState(), local_deferredLogger);
547 for (
auto& well : well_container_) {
548 if (well->isVFPActive(local_deferredLogger)){
549 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
555 this->updateWellPotentials(reportStepIdx,
557 simulator_.vanguard().summaryConfig(),
558 local_deferredLogger);
559 }
catch ( std::runtime_error& e ) {
560 const std::string msg =
"A zero well potential is returned for output purposes. ";
561 local_deferredLogger.warning(
"WELL_POTENTIAL_CALCULATION_FAILED", msg);
565 const auto& comm = simulator_.vanguard().grid().comm();
566 std::vector<Scalar> pot(this->numPhases(), 0.0);
567 const Group& fieldGroup = this->schedule().getGroup(
"FIELD", reportStepIdx);
568 WellGroupHelpers<Scalar>::updateGuideRates(fieldGroup,
579 local_deferredLogger);
581 auto exc_type = ExceptionType::NONE;
583 if (this->schedule_[reportStepIdx].has_gpmaint()) {
584 for (
auto& calculator : regionalAveragePressureCalculator_) {
585 calculator.second->template defineState<ElementContext>(simulator_);
587 const double dt = simulator_.timeStepSize();
588 WellGroupHelpers<Scalar>::updateGpMaintTargetForGroups(fieldGroup,
590 regionalAveragePressureCalculator_,
598 for (
auto& well : well_container_) {
599 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
600 + ScheduleEvents::INJECTION_TYPE_CHANGED
601 + ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
602 + ScheduleEvents::NEW_WELL;
604 const auto& events = this->schedule()[reportStepIdx].wellgroup_events();
605 const bool event = this->report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
606 const bool dyn_status_change = this->wellState().well(well->name()).status
607 != this->prevWellState().well(well->name()).status;
609 if (event || dyn_status_change) {
611 well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), local_deferredLogger);
612 well->calculateExplicitQuantities(simulator_, this->wellState(), local_deferredLogger);
613 well->solveWellEquation(simulator_, this->wellState(), this->groupState(), local_deferredLogger);
614 }
catch (
const std::exception& e) {
615 const std::string msg =
"Compute initial well solution for new well " + well->name() +
" failed. Continue with zero initial rates";
616 local_deferredLogger.warning(
"WELL_INITIAL_SOLVE_FAILED", msg);
622 OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
624 if (exc_type != ExceptionType::NONE) {
625 const std::string msg =
"Compute initial well solution for new wells failed. Continue with zero initial rates";
626 local_deferredLogger.warning(
"WELL_INITIAL_SOLVE_FAILED", msg);
629 logAndCheckForExceptionsAndThrow(local_deferredLogger,
630 exc_type,
"beginTimeStep() failed: " + exc_msg, this->terminal_output_, comm);
634 template<
typename TypeTag>
636 BlackoilWellModel<TypeTag>::wellTesting(
const int timeStepIdx,
637 const double simulationTime,
638 DeferredLogger& deferred_logger)
640 for (
const std::string& well_name : this->getWellsForTesting(timeStepIdx, simulationTime)) {
641 const Well& wellEcl = this->schedule().getWell(well_name, timeStepIdx);
642 if (wellEcl.getStatus() == Well::Status::SHUT)
645 WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
647 well->init(&this->phase_usage_, depth_, gravity_, B_avg_,
true);
649 Scalar well_efficiency_factor = wellEcl.getEfficiencyFactor();
650 WellGroupHelpers<Scalar>::accumulateGroupEfficiencyFactor(this->schedule().getGroup(wellEcl.groupName(),
654 well_efficiency_factor);
656 well->setWellEfficiencyFactor(well_efficiency_factor);
657 well->setVFPProperties(this->vfp_properties_.get());
658 well->setGuideRate(&this->guideRate_);
661 if (well->isProducer()) {
662 well->updateWellStateRates(simulator_, this->wellState(), deferred_logger);
664 if (well->isVFPActive(deferred_logger)) {
665 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
669 well->wellTesting(simulator_, simulationTime, this->wellState(),
670 this->groupState(), this->wellTestState(), deferred_logger);
671 }
catch (
const std::exception& e) {
672 const std::string msg = fmt::format(
"Exception during testing of well: {}. The well will not open.\n Exception message: {}", wellEcl.name(), e.what());
673 deferred_logger.warning(
"WELL_TESTING_FAILED", msg);
683 template<
typename TypeTag>
685 BlackoilWellModel<TypeTag>::
689 for (
auto&& pinfo : this->local_parallel_well_info_)
700 template<
typename TypeTag>
701 const SimulatorReportSingle&
702 BlackoilWellModel<TypeTag>::
703 lastReport()
const {
return last_report_; }
710 template<
typename TypeTag>
712 BlackoilWellModel<TypeTag>::
713 timeStepSucceeded(
const double simulationTime,
const double dt)
715 this->closed_this_step_.clear();
718 this->report_step_starts_ =
false;
719 const int reportStepIdx = simulator_.episodeIndex();
721 DeferredLogger local_deferredLogger;
722 for (
const auto& well : well_container_) {
723 if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
724 well->updateWaterThroughput(dt, this->wellState());
728 for (
const auto& well : well_container_) {
729 well->updateConnectionTransmissibilityFactor(simulator_, this->wellState().well(well->indexOfWell()));
730 well->updateConnectionDFactor(simulator_, this->wellState().well(well->indexOfWell()));
733 if (Indices::waterEnabled) {
734 this->updateFiltrationModelsPostStep(dt, FluidSystem::waterPhaseIdx, local_deferredLogger);
738 this->updateInjMult(local_deferredLogger);
741 for (
const auto& well : well_container_) {
742 well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);
745 if (this->terminal_output_) {
747 for (
const auto& [name, to] : this->switched_prod_groups_) {
748 const Group::ProductionCMode& oldControl = this->prevWGState().group_state.production_control(name);
749 std::string from = Group::ProductionCMode2String(oldControl);
751 std::string msg =
" Production Group " + name
752 +
" control mode changed from ";
755 local_deferredLogger.info(msg);
758 for (
const auto& [key, to] : this->switched_inj_groups_) {
759 const std::string& name = key.first;
760 const Opm::Phase& phase = key.second;
762 const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(name, phase);
763 std::string from = Group::InjectionCMode2String(oldControl);
765 std::string msg =
" Injection Group " + name
766 +
" control mode changed from ";
769 local_deferredLogger.info(msg);
775 rateConverter_->template defineState<ElementContext>(simulator_);
779 this->updateWellPotentials(reportStepIdx,
781 simulator_.vanguard().summaryConfig(),
782 local_deferredLogger);
783 }
catch ( std::runtime_error& e ) {
784 const std::string msg =
"A zero well potential is returned for output purposes. ";
785 local_deferredLogger.warning(
"WELL_POTENTIAL_CALCULATION_FAILED", msg);
788 updateWellTestState(simulationTime, this->wellTestState());
791 const Group& fieldGroup = this->schedule_.getGroup(
"FIELD", reportStepIdx);
792 this->checkGEconLimits(fieldGroup, simulationTime,
793 simulator_.episodeIndex(), local_deferredLogger);
794 this->checkGconsaleLimits(fieldGroup, this->wellState(),
795 simulator_.episodeIndex(), local_deferredLogger);
797 this->calculateProductivityIndexValues(local_deferredLogger);
799 this->commitWGState();
801 const Opm::Parallel::Communication& comm = grid().comm();
803 if (this->terminal_output_) {
804 global_deferredLogger.logMessages();
808 this->computeWellTemperature();
812 template<
typename TypeTag>
814 BlackoilWellModel<TypeTag>::
815 computeTotalRatesForDof(RateVector& rate,
816 unsigned elemIdx)
const
820 if (!is_cell_perforated_[elemIdx])
823 for (
const auto& well : well_container_)
824 well->addCellRates(rate, elemIdx);
828 template<
typename TypeTag>
829 template <
class Context>
831 BlackoilWellModel<TypeTag>::
832 computeTotalRatesForDof(RateVector& rate,
833 const Context& context,
835 unsigned timeIdx)
const
838 int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
840 if (!is_cell_perforated_[elemIdx])
843 for (
const auto& well : well_container_)
844 well->addCellRates(rate, elemIdx);
849 template<
typename TypeTag>
851 BlackoilWellModel<TypeTag>::
852 initializeWellState(
const int timeStepIdx)
854 const auto pressIx = []()
856 if (Indices::oilEnabled) {
return FluidSystem::oilPhaseIdx; }
857 if (Indices::waterEnabled) {
return FluidSystem::waterPhaseIdx; }
859 return FluidSystem::gasPhaseIdx;
862 auto cellPressures = std::vector<Scalar>(this->local_num_cells_, Scalar{0});
863 auto cellTemperatures = std::vector<Scalar>(this->local_num_cells_, Scalar{0});
865 auto elemCtx = ElementContext { this->simulator_ };
866 const auto& gridView = this->simulator_.vanguard().gridView();
868 OPM_BEGIN_PARALLEL_TRY_CATCH();
869 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
870 elemCtx.updatePrimaryStencil(elem);
871 elemCtx.updatePrimaryIntensiveQuantities(0);
873 const auto ix = elemCtx.globalSpaceIndex(0, 0);
874 const auto& fs = elemCtx.intensiveQuantities(0, 0).fluidState();
876 cellPressures[ix] = fs.pressure(pressIx).value();
877 cellTemperatures[ix] = fs.temperature(0).value();
879 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel::initializeWellState() failed: ",
880 this->simulator_.vanguard().grid().comm());
882 this->wellState().init(cellPressures, cellTemperatures, this->schedule(), this->wells_ecl_,
883 this->local_parallel_well_info_, timeStepIdx,
884 &this->prevWellState(), this->well_perf_data_,
885 this->summaryState());
892 template<
typename TypeTag>
894 BlackoilWellModel<TypeTag>::
895 createWellContainer(
const int report_step)
897 DeferredLogger local_deferredLogger;
899 const int nw = this->numLocalWells();
901 well_container_.clear();
904 well_container_.reserve(nw);
906 for (
int w = 0; w < nw; ++w) {
907 const Well& well_ecl = this->wells_ecl_[w];
909 if (!well_ecl.hasConnections()) {
914 const std::string& well_name = well_ecl.name();
915 const auto well_status = this->schedule()
916 .getWell(well_name, report_step).getStatus();
918 if ((well_ecl.getStatus() == Well::Status::SHUT) ||
919 (well_status == Well::Status::SHUT))
922 if (well_ecl.getStatus() != Well::Status::SHUT) {
923 this->closed_this_step_.insert(well_name);
924 this->wellState().shutWell(w);
931 if (this->wellTestState().well_is_closed(well_name)) {
936 const bool closed_this_step = (this->wellTestState().lastTestTime(well_name) == simulator_.time());
939 auto& events = this->wellState().well(w).events;
940 if (events.hasEvent(ScheduleEvents::REQUEST_OPEN_WELL)) {
941 if (!closed_this_step) {
942 this->wellTestState().open_well(well_name);
943 this->wellTestState().open_completions(well_name);
945 events.clearEvent(ScheduleEvents::REQUEST_OPEN_WELL);
951 bool wellIsStopped =
false;
952 if (this->wellTestState().well_is_closed(well_name))
954 if (well_ecl.getAutomaticShutIn()) {
956 this->wellState().shutWell(w);
959 if (!well_ecl.getAllowCrossFlow()) {
962 this->wellState().shutWell(w);
966 this->wellState().stopWell(w);
967 wellIsStopped =
true;
972 if (!well_ecl.getAllowCrossFlow()) {
973 const bool any_zero_rate_constraint = well_ecl.isProducer()
974 ? well_ecl.productionControls(this->summaryState_).anyZeroRateConstraint()
975 : well_ecl.injectionControls(this->summaryState_).anyZeroRateConstraint();
976 if (any_zero_rate_constraint) {
978 local_deferredLogger.debug(fmt::format(
" Well {} gets shut due to having zero rate constraint and disallowing crossflow ", well_ecl.name()) );
979 this->wellState().shutWell(w);
984 if (well_status == Well::Status::STOP) {
985 this->wellState().stopWell(w);
986 wellIsStopped =
true;
989 well_container_.emplace_back(this->createWellPointer(w, report_step));
992 well_container_.back()->stopWell();
998 const Opm::Parallel::Communication& comm = grid().comm();
1000 if (this->terminal_output_) {
1001 global_deferredLogger.logMessages();
1004 this->well_container_generic_.clear();
1005 for (
auto& w : well_container_)
1006 this->well_container_generic_.push_back(w.get());
1008 const auto& network = this->schedule()[report_step].network();
1009 if (network.active() && !this->node_pressures_.empty()) {
1010 for (
auto& well: this->well_container_generic_) {
1014 if (well->isProducer()) {
1015 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
1016 if (it != this->node_pressures_.end()) {
1019 const Scalar nodal_pressure = it->second;
1020 well->setDynamicThpLimit(nodal_pressure);
1026 this->registerOpenWellsForWBPCalculation();
1033 template <
typename TypeTag>
1034 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1035 BlackoilWellModel<TypeTag>::
1036 createWellPointer(
const int wellID,
const int report_step)
const
1038 const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
1040 if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
1041 return this->
template createTypedWellPointer<StandardWell<TypeTag>>(wellID, report_step);
1044 return this->
template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, report_step);
1052 template <
typename TypeTag>
1053 template <
typename WellType>
1054 std::unique_ptr<WellType>
1055 BlackoilWellModel<TypeTag>::
1056 createTypedWellPointer(
const int wellID,
const int time_step)
const
1059 const auto& perf_data = this->well_perf_data_[wellID];
1062 const auto pvtreg = perf_data.empty()
1063 ? 0 : this->pvt_region_idx_[perf_data.front().cell_index];
1065 const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get();
1066 const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
1068 return std::make_unique<WellType>(this->wells_ecl_[wellID],
1072 *this->rateConverter_,
1074 this->numComponents(),
1084 template<
typename TypeTag>
1085 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1086 BlackoilWellModel<TypeTag>::
1087 createWellForWellTest(
const std::string& well_name,
1088 const int report_step,
1089 DeferredLogger& deferred_logger)
const
1092 const int nw_wells_ecl = this->wells_ecl_.size();
1093 int index_well_ecl = 0;
1094 for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
1095 if (well_name == this->wells_ecl_[index_well_ecl].name()) {
1100 if (index_well_ecl == nw_wells_ecl) {
1101 OPM_DEFLOG_THROW(std::logic_error,
1102 fmt::format(
"Could not find well {} in wells_ecl ", well_name),
1106 return this->createWellPointer(index_well_ecl, report_step);
1111 template<
typename TypeTag>
1113 BlackoilWellModel<TypeTag>::
1114 doPreStepNetworkRebalance(DeferredLogger& deferred_logger) {
1115 const double dt = this->simulator_.timeStepSize();
1117 auto& well_state = this->wellState();
1118 const std::size_t max_iter = param_.network_max_iterations_;
1119 bool converged =
false;
1120 std::size_t iter = 0;
1121 bool changed_well_group =
false;
1123 changed_well_group = updateWellControlsAndNetwork(
true, dt, deferred_logger);
1124 assembleWellEqWithoutIteration(dt, deferred_logger);
1125 converged = this->getWellConvergence(this->B_avg_,
true).converged() && !changed_well_group;
1130 for (
auto& well : this->well_container_) {
1131 well->solveEqAndUpdateWellState(simulator_, well_state, deferred_logger);
1133 this->initPrimaryVariablesEvaluation();
1134 }
while (iter < max_iter);
1137 const std::string msg = fmt::format(
"Initial (pre-step) network balance did not get converged with {} iterations, "
1138 "unconverged network balance result will be used", max_iter);
1139 deferred_logger.warning(msg);
1141 const std::string msg = fmt::format(
"Initial (pre-step) network balance converged with {} iterations", iter);
1142 deferred_logger.debug(msg);
1149 template<
typename TypeTag>
1151 BlackoilWellModel<TypeTag>::
1152 assemble(
const int iterationIdx,
1156 DeferredLogger local_deferredLogger;
1157 if (this->glift_debug) {
1158 const std::string msg = fmt::format(
1159 "assemble() : iteration {}" , iterationIdx);
1160 this->gliftDebug(msg, local_deferredLogger);
1162 last_report_ = SimulatorReportSingle();
1163 Dune::Timer perfTimer;
1165 this->closed_offending_wells_.clear();
1168 const int episodeIdx = simulator_.episodeIndex();
1169 const auto& network = this->schedule()[episodeIdx].network();
1170 if (!this->wellsActive() && !network.active()) {
1175 if (iterationIdx == 0 && this->wellsActive()) {
1180 OPM_BEGIN_PARALLEL_TRY_CATCH();
1182 calculateExplicitQuantities(local_deferredLogger);
1183 prepareTimeStep(local_deferredLogger);
1185 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1186 "assemble() failed (It=0): ",
1187 this->terminal_output_, grid().comm());
1190 const bool well_group_control_changed = updateWellControlsAndNetwork(
false, dt, local_deferredLogger);
1194 if ( ! this->wellsActive() ) {
1198 assembleWellEqWithoutIteration(dt, local_deferredLogger);
1202 last_report_.well_group_control_changed = well_group_control_changed;
1203 last_report_.assemble_time_well += perfTimer.stop();
1209 template<
typename TypeTag>
1211 BlackoilWellModel<TypeTag>::
1212 updateWellControlsAndNetwork(
const bool mandatory_network_balance,
const double dt, DeferredLogger& local_deferredLogger)
1215 bool do_network_update =
true;
1216 bool well_group_control_changed =
false;
1218 const std::size_t iteration_to_relax = param_.network_max_strict_iterations_;
1220 const std::size_t max_iteration = param_.network_max_iterations_;
1221 std::size_t network_update_iteration = 0;
1222 while (do_network_update) {
1223 if (this->terminal_output_ && (network_update_iteration == iteration_to_relax) ) {
1224 local_deferredLogger.info(
" we begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) +
" iterations ");
1226 const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
1227 std::tie(do_network_update, well_group_control_changed) =
1228 updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger);
1229 ++network_update_iteration;
1231 if (network_update_iteration >= max_iteration ) {
1232 if (this->terminal_output_) {
1233 local_deferredLogger.info(
"maximum of " + std::to_string(max_iteration) +
" iterations has been used, we stop the network update now. "
1234 "The simulation will continue with unconverged network results");
1239 return well_group_control_changed;
1245 template<
typename TypeTag>
1246 std::pair<bool, bool>
1247 BlackoilWellModel<TypeTag>::
1248 updateWellControlsAndNetworkIteration(
const bool mandatory_network_balance,
1249 const bool relax_network_tolerance,
1251 DeferredLogger& local_deferredLogger)
1253 auto [well_group_control_changed, more_network_update] =
1254 updateWellControls(mandatory_network_balance, local_deferredLogger, relax_network_tolerance);
1256 bool alq_updated =
false;
1257 OPM_BEGIN_PARALLEL_TRY_CATCH();
1260 initPrimaryVariablesEvaluation();
1262 alq_updated = maybeDoGasLiftOptimize(local_deferredLogger);
1264 prepareWellsBeforeAssembling(dt, local_deferredLogger);
1266 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
"updateWellControlsAndNetworkIteration() failed: ",
1267 this->terminal_output_, grid().comm());
1270 const int reportStepIdx = simulator_.episodeIndex();
1271 if (alq_updated || BlackoilWellModelGuideRates(*this).
1272 guideRateUpdateIsNeeded(reportStepIdx)) {
1273 const double simulationTime = simulator_.time();
1274 const auto& comm = simulator_.vanguard().grid().comm();
1275 const auto& summaryState = simulator_.vanguard().summaryState();
1276 std::vector<Scalar> pot(this->numPhases(), 0.0);
1277 const Group& fieldGroup = this->schedule().getGroup(
"FIELD", reportStepIdx);
1278 WellGroupHelpers<Scalar>::updateGuideRates(fieldGroup,
1289 local_deferredLogger);
1292 return {more_network_update, well_group_control_changed};
1298 template <
typename TypeTag>
1300 BlackoilWellModel<TypeTag>::
1301 computeWellGroupThp(
const double dt, DeferredLogger& local_deferredLogger)
1303 const int reportStepIdx = this->simulator_.episodeIndex();
1304 const auto& network = this->schedule()[reportStepIdx].network();
1305 const auto& balance = this->schedule()[reportStepIdx].network_balance();
1306 const Scalar thp_tolerance = balance.thp_tolerance();
1308 if (!network.active()) {
1312 auto& well_state = this->wellState();
1313 auto& group_state = this->groupState();
1315 for (
const std::string& nodeName : network.node_names()) {
1316 const bool has_choke = network.node(nodeName).as_choke();
1318 const auto& summary_state = this->simulator_.vanguard().summaryState();
1319 const Group& group = this->schedule().getGroup(nodeName, reportStepIdx);
1320 const auto ctrl = group.productionControls(summary_state);
1321 const auto cmode = ctrl.cmode;
1322 const auto pu = this->phase_usage_;
1324 std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
1325 Scalar gratTargetFromSales = 0.0;
1326 if (group_state.has_grat_sales_target(group.name()))
1327 gratTargetFromSales = group_state.grat_sales_target(group.name());
1329 WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
1330 gratTargetFromSales, nodeName, group_state,
1331 group.has_gpmaint_control(cmode));
1332 const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
1334 auto mismatch = [&] (
auto group_thp) {
1335 Scalar group_rate(0.0);
1337 for (
auto& well : this->well_container_) {
1338 std::string well_name = well->name();
1339 auto& ws = well_state.well(well_name);
1340 if (group.hasWell(well_name)) {
1341 well->setDynamicThpLimit(group_thp);
1342 const Well& well_ecl = this->wells_ecl_[well->indexOfWell()];
1343 const auto inj_controls = Well::InjectionControls(0);
1344 const auto prod_controls = well_ecl.productionControls(summary_state);
1345 well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger,
false,
false);
1346 rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
1350 return (group_rate - orig_target)/orig_target;
1353 const auto upbranch = network.uptree_branch(nodeName);
1354 const auto it = this->node_pressures_.find((*upbranch).uptree_node());
1355 const Scalar nodal_pressure = it->second;
1356 Scalar well_group_thp = nodal_pressure;
1358 std::optional<Scalar> autochoke_thp;
1359 if (
auto iter = this->well_group_thp_calc_.find(nodeName); iter != this->well_group_thp_calc_.end()) {
1360 autochoke_thp = this->well_group_thp_calc_.at(nodeName);
1364 std::array<Scalar, 2> range_initial;
1365 if (!autochoke_thp.has_value()){
1366 Scalar min_thp, max_thp;
1368 std::string node_name = nodeName;
1369 while (!network.node(node_name).terminal_pressure().has_value()) {
1370 auto branch = network.uptree_branch(node_name).value();
1371 node_name = branch.uptree_node();
1373 min_thp = network.node(node_name).terminal_pressure().value();
1374 WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp);
1377 std::array<Scalar, 2> range = {Scalar{0.9}*min_thp, Scalar{1.1}*max_thp};
1378 std::optional<Scalar> appr_sol;
1379 WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger);
1382 range_initial = {min_thp, max_thp};
1385 if (!autochoke_thp.has_value() || autochoke_thp.value() > nodal_pressure) {
1387 std::array<Scalar, 2> range = autochoke_thp.has_value() ?
1388 std::array<Scalar, 2>{Scalar{0.9} * autochoke_thp.value(),
1389 Scalar{1.1} * autochoke_thp.value()} : range_initial;
1391 std::optional<Scalar> approximate_solution;
1392 const Scalar tolerance1 = thp_tolerance;
1393 local_deferredLogger.debug(
"Using brute force search to bracket the group THP");
1394 const bool finding_bracket = WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger);
1396 if (approximate_solution.has_value()) {
1397 autochoke_thp = *approximate_solution;
1398 local_deferredLogger.debug(
"Approximate group THP value found: " + std::to_string(autochoke_thp.value()));
1399 }
else if (finding_bracket) {
1400 const Scalar tolerance2 = thp_tolerance;
1401 const int max_iteration_solve = 100;
1403 autochoke_thp = RegulaFalsiBisection<ThrowOnError>::
1404 solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration);
1405 local_deferredLogger.debug(
" bracket = [" + std::to_string(low) +
", " + std::to_string(high) +
"], " +
1406 "iteration = " + std::to_string(iteration));
1407 local_deferredLogger.debug(
"Group THP value = " + std::to_string(autochoke_thp.value()));
1409 autochoke_thp.reset();
1410 local_deferredLogger.debug(
"Group THP solve failed due to bracketing failure");
1413 if (autochoke_thp.has_value()) {
1414 well_group_thp_calc_[nodeName] = autochoke_thp.value();
1417 well_group_thp = std::max(autochoke_thp.value(), nodal_pressure);
1420 for (
auto& well : this->well_container_) {
1421 std::string well_name = well->name();
1422 if (group.hasWell(well_name)) {
1423 well->setDynamicThpLimit(well_group_thp);
1428 group_state.update_well_group_thp(nodeName, well_group_thp);
1433 template<
typename TypeTag>
1435 BlackoilWellModel<TypeTag>::
1436 assembleDomain([[maybe_unused]]
const int iterationIdx,
1438 const Domain& domain)
1440 last_report_ = SimulatorReportSingle();
1441 Dune::Timer perfTimer;
1445 const int episodeIdx = simulator_.episodeIndex();
1446 const auto& network = this->schedule()[episodeIdx].network();
1447 if (!this->wellsActive() && !network.active()) {
1457 DeferredLogger local_deferredLogger;
1458 updateWellControlsDomain(local_deferredLogger, domain);
1459 initPrimaryVariablesEvaluationDomain(domain);
1460 assembleWellEqDomain(dt, domain, local_deferredLogger);
1464 if (this->terminal_output_) {
1465 local_deferredLogger.logMessages();
1468 last_report_.converged =
true;
1469 last_report_.assemble_time_well += perfTimer.stop();
1473 template<
typename TypeTag>
1475 BlackoilWellModel<TypeTag>::
1476 maybeDoGasLiftOptimize(DeferredLogger& deferred_logger)
1478 bool do_glift_optimization =
false;
1479 int num_wells_changed = 0;
1480 const double simulation_time = simulator_.time();
1481 const Scalar min_wait = simulator_.vanguard().schedule().glo(simulator_.episodeIndex()).min_wait();
1486 if ( simulation_time == this->last_glift_opt_time_ || simulation_time >= (this->last_glift_opt_time_ + min_wait)) {
1487 do_glift_optimization =
true;
1488 this->last_glift_opt_time_ = simulation_time;
1491 if (do_glift_optimization) {
1492 GLiftOptWells glift_wells;
1493 GLiftProdWells prod_wells;
1494 GLiftWellStateMap state_map;
1502 GLiftEclWells ecl_well_map;
1503 initGliftEclWellMap(ecl_well_map);
1504 GasLiftGroupInfo group_info {
1506 simulator_.vanguard().schedule(),
1507 simulator_.vanguard().summaryState(),
1508 simulator_.episodeIndex(),
1509 simulator_.model().newtonMethod().numIterations(),
1514 simulator_.vanguard().grid().comm(),
1517 group_info.initialize();
1518 gasLiftOptimizationStage1(deferred_logger, prod_wells, glift_wells,
1519 group_info, state_map);
1520 this->gasLiftOptimizationStage2(deferred_logger, prod_wells, glift_wells,
1521 group_info, state_map, simulator_.episodeIndex());
1522 if (this->glift_debug) {
1523 this->gliftDebugShowALQ(deferred_logger);
1525 num_wells_changed = glift_wells.size();
1527 num_wells_changed = this->comm_.sum(num_wells_changed);
1528 return num_wells_changed > 0;
1531 template<
typename TypeTag>
1533 BlackoilWellModel<TypeTag>::
1534 gasLiftOptimizationStage1(DeferredLogger& deferred_logger,
1535 GLiftProdWells& prod_wells,
1536 GLiftOptWells &glift_wells,
1537 GasLiftGroupInfo<Scalar>& group_info,
1538 GLiftWellStateMap& state_map)
1540 auto comm = simulator_.vanguard().grid().comm();
1541 int num_procs = comm.size();
1567 for (
int i = 0; i< num_procs; i++) {
1568 int num_rates_to_sync = 0;
1569 GLiftSyncGroups groups_to_sync;
1570 if (comm.rank() == i) {
1572 for (
const auto& well : well_container_) {
1574 if (group_info.hasWell(well->name())) {
1575 gasLiftOptimizationStage1SingleWell(
1576 well.get(), deferred_logger, prod_wells, glift_wells,
1577 group_info, state_map, groups_to_sync
1581 num_rates_to_sync = groups_to_sync.size();
1583 num_rates_to_sync = comm.sum(num_rates_to_sync);
1584 if (num_rates_to_sync > 0) {
1585 std::vector<int> group_indexes;
1586 group_indexes.reserve(num_rates_to_sync);
1587 std::vector<Scalar> group_alq_rates;
1588 group_alq_rates.reserve(num_rates_to_sync);
1589 std::vector<Scalar> group_oil_rates;
1590 group_oil_rates.reserve(num_rates_to_sync);
1591 std::vector<Scalar> group_gas_rates;
1592 group_gas_rates.reserve(num_rates_to_sync);
1593 std::vector<Scalar> group_water_rates;
1594 group_water_rates.reserve(num_rates_to_sync);
1595 if (comm.rank() == i) {
1596 for (
auto idx : groups_to_sync) {
1597 auto [oil_rate, gas_rate, water_rate, alq] = group_info.getRates(idx);
1598 group_indexes.push_back(idx);
1599 group_oil_rates.push_back(oil_rate);
1600 group_gas_rates.push_back(gas_rate);
1601 group_water_rates.push_back(water_rate);
1602 group_alq_rates.push_back(alq);
1605 group_indexes.resize(num_rates_to_sync);
1606 group_oil_rates.resize(num_rates_to_sync);
1607 group_gas_rates.resize(num_rates_to_sync);
1608 group_water_rates.resize(num_rates_to_sync);
1609 group_alq_rates.resize(num_rates_to_sync);
1612 Parallel::MpiSerializer ser(comm);
1613 ser.broadcast(i, group_indexes, group_oil_rates,
1614 group_gas_rates, group_water_rates, group_alq_rates);
1616 if (comm.rank() != i) {
1617 for (
int j=0; j<num_rates_to_sync; j++) {
1618 group_info.updateRate(group_indexes[j],
1619 group_oil_rates[j], group_gas_rates[j], group_water_rates[j], group_alq_rates[j]);
1622 if (this->glift_debug) {
1624 if (comm.rank() == i) {
1625 counter = this->wellState().gliftGetDebugCounter();
1627 counter = comm.sum(counter);
1628 if (comm.rank() != i) {
1629 this->wellState().gliftSetDebugCounter(counter);
1639 template<
typename TypeTag>
1641 BlackoilWellModel<TypeTag>::
1642 gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag>* well,
1643 DeferredLogger& deferred_logger,
1644 GLiftProdWells& prod_wells,
1645 GLiftOptWells& glift_wells,
1646 GasLiftGroupInfo<Scalar>& group_info,
1647 GLiftWellStateMap& state_map,
1648 GLiftSyncGroups& sync_groups)
1650 const auto& summary_state = simulator_.vanguard().summaryState();
1651 std::unique_ptr<GasLiftSingleWell> glift
1652 = std::make_unique<GasLiftSingleWell>(
1653 *well, simulator_, summary_state,
1654 deferred_logger, this->wellState(), this->groupState(),
1655 group_info, sync_groups, this->comm_, this->glift_debug);
1656 auto state = glift->runOptimize(
1657 simulator_.model().newtonMethod().numIterations());
1659 state_map.insert({well->name(), std::move(state)});
1660 glift_wells.insert({well->name(), std::move(glift)});
1663 prod_wells.insert({well->name(), well});
1667 template<
typename TypeTag>
1669 BlackoilWellModel<TypeTag>::
1670 initGliftEclWellMap(GLiftEclWells &ecl_well_map)
1672 for (
const auto& well: well_container_ ) {
1673 ecl_well_map.try_emplace(
1674 well->name(), &(well->wellEcl()), well->indexOfWell());
1679 template<
typename TypeTag>
1681 BlackoilWellModel<TypeTag>::
1682 assembleWellEq(
const double dt, DeferredLogger& deferred_logger)
1684 for (
auto& well : well_container_) {
1685 well->assembleWellEq(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1690 template<
typename TypeTag>
1692 BlackoilWellModel<TypeTag>::
1693 assembleWellEqDomain(
const double dt,
const Domain& domain, DeferredLogger& deferred_logger)
1695 for (
auto& well : well_container_) {
1696 if (well_domain_.at(well->name()) == domain.index) {
1697 well->assembleWellEq(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1703 template<
typename TypeTag>
1705 BlackoilWellModel<TypeTag>::
1706 prepareWellsBeforeAssembling(
const double dt, DeferredLogger& deferred_logger)
1708 for (
auto& well : well_container_) {
1709 well->prepareWellBeforeAssembling(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1714 template<
typename TypeTag>
1716 BlackoilWellModel<TypeTag>::
1717 assembleWellEqWithoutIteration(
const double dt, DeferredLogger& deferred_logger)
1721 OPM_BEGIN_PARALLEL_TRY_CATCH();
1723 for (
auto& well: well_container_) {
1724 well->assembleWellEqWithoutIteration(simulator_, dt, this->wellState(), this->groupState(),
1727 OPM_END_PARALLEL_TRY_CATCH_LOG(deferred_logger,
"BlackoilWellModel::assembleWellEqWithoutIteration failed: ",
1728 this->terminal_output_, grid().comm());
1733 template<
typename TypeTag>
1735 BlackoilWellModel<TypeTag>::
1736 apply(
const BVector& x, BVector& Ax)
const
1738 for (
auto& well : well_container_) {
1740 const auto& cells = well->cells();
1741 x_local_.resize(cells.size());
1742 Ax_local_.resize(cells.size());
1744 for (
size_t i = 0; i < cells.size(); ++i) {
1745 x_local_[i] = x[cells[i]];
1746 Ax_local_[i] = Ax[cells[i]];
1749 well->apply(x_local_, Ax_local_);
1751 for (
size_t i = 0; i < cells.size(); ++i) {
1753 Ax[cells[i]] = Ax_local_[i];
1759 template<
typename TypeTag>
1761 BlackoilWellModel<TypeTag>::
1762 applyDomain(
const BVector& x, BVector& Ax,
const int domainIndex)
const
1764 for (
size_t well_index = 0; well_index < well_container_.size(); ++well_index) {
1765 auto& well = well_container_[well_index];
1766 if (well_domain_.at(well->name()) == domainIndex) {
1769 const auto& local_cells = well_local_cells_[well_index];
1770 x_local_.resize(local_cells.size());
1771 Ax_local_.resize(local_cells.size());
1773 for (
size_t i = 0; i < local_cells.size(); ++i) {
1774 x_local_[i] = x[local_cells[i]];
1775 Ax_local_[i] = Ax[local_cells[i]];
1778 well->apply(x_local_, Ax_local_);
1780 for (
size_t i = 0; i < local_cells.size(); ++i) {
1782 Ax[local_cells[i]] = Ax_local_[i];
1788#if COMPILE_BDA_BRIDGE
1789 template<
typename TypeTag>
1791 BlackoilWellModel<TypeTag>::
1792 getWellContributions(WellContributions<Scalar>& wellContribs)
const
1795 wellContribs.setBlockSize(StandardWell<TypeTag>::Indices::numEq, StandardWell<TypeTag>::numStaticWellEq);
1797 for(
unsigned int i = 0; i < well_container_.size(); i++){
1798 auto& well = well_container_[i];
1799 std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1801 wellContribs.addNumBlocks(derived->linSys().getNumBlocks());
1806 wellContribs.alloc();
1808 for(
unsigned int i = 0; i < well_container_.size(); i++){
1809 auto& well = well_container_[i];
1811 auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag>>(well);
1813 derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs);
1815 auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
1817 derived_ms->linSys().extract(wellContribs);
1819 OpmLog::warning(
"Warning unknown type of well");
1827 template<
typename TypeTag>
1829 BlackoilWellModel<TypeTag>::
1830 applyScaleAdd(
const Scalar alpha,
const BVector& x, BVector& Ax)
const
1832 if (this->well_container_.empty()) {
1836 if( scaleAddRes_.size() != Ax.size() ) {
1837 scaleAddRes_.resize( Ax.size() );
1842 apply( x, scaleAddRes_ );
1844 Ax.axpy( alpha, scaleAddRes_ );
1848 template<
typename TypeTag>
1850 BlackoilWellModel<TypeTag>::
1851 applyScaleAddDomain(
const Scalar alpha,
const BVector& x, BVector& Ax,
const int domainIndex)
const
1853 if (this->well_container_.empty()) {
1857 if( scaleAddRes_.size() != Ax.size() ) {
1858 scaleAddRes_.resize( Ax.size() );
1863 applyDomain(x, scaleAddRes_, domainIndex);
1865 Ax.axpy( alpha, scaleAddRes_ );
1868 template<
typename TypeTag>
1870 BlackoilWellModel<TypeTag>::
1871 addWellContributions(SparseMatrixAdapter& jacobian)
const
1873 for (
const auto& well: well_container_ ) {
1874 well->addWellContributions(jacobian);
1878 template<
typename TypeTag>
1880 BlackoilWellModel<TypeTag>::
1881 addWellPressureEquations(PressureMatrix& jacobian,
const BVector& weights,
const bool use_well_weights)
const
1883 int nw = this->numLocalWellsEnd();
1884 int rdofs = local_num_cells_;
1885 for (
int i = 0; i < nw; i++ ) {
1886 int wdof = rdofs + i;
1887 jacobian[wdof][wdof] = 1.0;
1890 for (
const auto& well : well_container_ ) {
1891 well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
1895 template<
typename TypeTag>
1897 BlackoilWellModel<TypeTag>::
1898 addWellPressureEquationsDomain([[maybe_unused]] PressureMatrix& jacobian,
1899 [[maybe_unused]]
const BVector& weights,
1900 [[maybe_unused]]
const bool use_well_weights,
1901 [[maybe_unused]]
const int domainIndex)
const
1903 throw std::logic_error(
"CPRW is not yet implemented for NLDD subdomains");
1920 template <
typename TypeTag>
1921 void BlackoilWellModel<TypeTag>::
1922 addReservoirSourceTerms(GlobalEqVector& residual,
1923 std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress)
const
1928 for (
const auto& well : well_container_) {
1929 if (!well->isOperableAndSolvable() && !well->wellIsStopped()) {
1932 const auto& cells = well->cells();
1933 const auto& rates = well->connectionRates();
1934 for (
unsigned perfIdx = 0; perfIdx < rates.size(); ++perfIdx) {
1935 unsigned cellIdx = cells[perfIdx];
1936 auto rate = rates[perfIdx];
1938 VectorBlockType res(0.0);
1939 using MatrixBlockType =
typename SparseMatrixAdapter::MatrixBlock;
1940 MatrixBlockType bMat(0.0);
1941 simulator_.model().linearizer().setResAndJacobi(res, bMat, rate);
1942 residual[cellIdx] += res;
1943 *diagMatAddress[cellIdx] += bMat;
1949 template<
typename TypeTag>
1951 BlackoilWellModel<TypeTag>::
1952 addWellPressureEquationsStruct(PressureMatrix& jacobian)
const
1954 int nw = this->numLocalWellsEnd();
1955 int rdofs = local_num_cells_;
1956 for(
int i=0; i < nw; i++){
1957 int wdof = rdofs + i;
1958 jacobian.entry(wdof,wdof) = 1.0;
1960 std::vector<std::vector<int>> wellconnections = this->getMaxWellConnections();
1961 for(
int i=0; i < nw; i++){
1962 const auto& perfcells = wellconnections[i];
1963 for(
int perfcell : perfcells){
1964 int wdof = rdofs + i;
1965 jacobian.entry(wdof,perfcell) = 0.0;
1966 jacobian.entry(perfcell, wdof) = 0.0;
1972 template<
typename TypeTag>
1974 BlackoilWellModel<TypeTag>::
1975 recoverWellSolutionAndUpdateWellState(
const BVector& x)
1977 DeferredLogger local_deferredLogger;
1978 OPM_BEGIN_PARALLEL_TRY_CATCH();
1980 for (
auto& well : well_container_) {
1981 const auto& cells = well->cells();
1982 x_local_.resize(cells.size());
1984 for (
size_t i = 0; i < cells.size(); ++i) {
1985 x_local_[i] = x[cells[i]];
1987 well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_, this->wellState(), local_deferredLogger);
1990 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1991 "recoverWellSolutionAndUpdateWellState() failed: ",
1992 this->terminal_output_, simulator_.vanguard().grid().comm());
1996 template<
typename TypeTag>
1998 BlackoilWellModel<TypeTag>::
1999 recoverWellSolutionAndUpdateWellStateDomain(
const BVector& x,
const Domain& domain)
2004 DeferredLogger local_deferredLogger;
2005 for (
auto& well : well_container_) {
2006 if (well_domain_.at(well->name()) == domain.index) {
2007 const auto& cells = well->cells();
2008 x_local_.resize(cells.size());
2010 for (
size_t i = 0; i < cells.size(); ++i) {
2011 x_local_[i] = x[cells[i]];
2013 well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_,
2015 local_deferredLogger);
2020 if (this->terminal_output_) {
2021 local_deferredLogger.logMessages();
2026 template<
typename TypeTag>
2028 BlackoilWellModel<TypeTag>::
2029 initPrimaryVariablesEvaluation()
const
2031 for (
auto& well : well_container_) {
2032 well->initPrimaryVariablesEvaluation();
2037 template<
typename TypeTag>
2039 BlackoilWellModel<TypeTag>::
2040 initPrimaryVariablesEvaluationDomain(
const Domain& domain)
const
2042 for (
auto& well : well_container_) {
2043 if (well_domain_.at(well->name()) == domain.index) {
2044 well->initPrimaryVariablesEvaluation();
2054 template<
typename TypeTag>
2056 BlackoilWellModel<TypeTag>::
2057 getDomainWellConvergence(
const Domain& domain,
2058 const std::vector<Scalar>& B_avg,
2059 DeferredLogger& local_deferredLogger)
const
2061 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
2062 const bool relax_tolerance = iterationIdx > param_.strict_outer_iter_wells_;
2064 ConvergenceReport report;
2065 for (
const auto& well : well_container_) {
2066 if ((well_domain_.at(well->name()) == domain.index)) {
2067 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
2068 report += well->getWellConvergence(simulator_,
2071 local_deferredLogger,
2074 ConvergenceReport xreport;
2075 using CR = ConvergenceReport;
2076 xreport.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
2083 if (this->terminal_output_) {
2084 for (
const auto& f : report.wellFailures()) {
2085 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
2086 local_deferredLogger.debug(
"NaN residual found with phase " + std::to_string(f.phase()) +
" for well " + f.wellName());
2087 }
else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
2088 local_deferredLogger.debug(
"Too large residual found with phase " + std::to_string(f.phase()) +
" for well " + f.wellName());
2099 template<
typename TypeTag>
2101 BlackoilWellModel<TypeTag>::
2102 getWellConvergence(
const std::vector<Scalar>& B_avg,
bool checkWellGroupControls)
const
2105 DeferredLogger local_deferredLogger;
2107 ConvergenceReport local_report;
2108 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
2109 for (
const auto& well : well_container_) {
2110 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
2111 local_report += well->getWellConvergence(
2112 simulator_, this->wellState(), B_avg, local_deferredLogger,
2113 iterationIdx > param_.strict_outer_iter_wells_);
2115 ConvergenceReport report;
2116 using CR = ConvergenceReport;
2117 report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
2118 local_report += report;
2122 const Opm::Parallel::Communication comm = grid().comm();
2127 if (checkWellGroupControls) {
2128 report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed);
2131 if (this->terminal_output_) {
2132 global_deferredLogger.logMessages();
2135 if (this->terminal_output_) {
2136 for (
const auto& f : report.wellFailures()) {
2137 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
2138 OpmLog::debug(
"NaN residual found with phase " + std::to_string(f.phase()) +
" for well " + f.wellName());
2139 }
else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
2140 OpmLog::debug(
"Too large residual found with phase " + std::to_string(f.phase()) +
" for well " + f.wellName());
2151 template<
typename TypeTag>
2157 for (
auto& well : well_container_) {
2158 well->calculateExplicitQuantities(simulator_, this->wellState(), deferred_logger);
2166 template<
typename TypeTag>
2167 std::pair<bool, bool>
2171 const int episodeIdx = simulator_.episodeIndex();
2172 const auto& network = this->schedule()[episodeIdx].network();
2173 if (!this->wellsActive() && !network.active()) {
2174 return {
false,
false};
2177 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
2178 const auto& comm = simulator_.vanguard().grid().comm();
2179 this->updateAndCommunicateGroupData(episodeIdx, iterationIdx);
2182 bool more_network_update =
false;
2183 if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
2184 const double dt = this->simulator_.timeStepSize();
2186 computeWellGroupThp(dt, deferred_logger);
2187 const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx);
2188 const Scalar network_imbalance = comm.max(local_network_imbalance);
2189 const auto& balance = this->schedule()[episodeIdx].network_balance();
2190 constexpr Scalar relaxation_factor = 10.0;
2191 const Scalar tolerance = relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() : balance.pressure_tolerance();
2192 more_network_update = this->networkActive() && network_imbalance > tolerance;
2195 bool changed_well_group =
false;
2197 const int nupcol = this->schedule()[episodeIdx].nupcol();
2200 if (iterationIdx <= nupcol) {
2201 const Group& fieldGroup = this->schedule().getGroup(
"FIELD", episodeIdx);
2202 changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
2205 bool changed_well_to_group =
false;
2209 OPM_BEGIN_PARALLEL_TRY_CATCH()
2210 for (const auto& well : well_container_) {
2211 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
2212 const bool changed_well = well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
2214 changed_well_to_group = changed_well || changed_well_to_group;
2217 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel: updating well controls failed: ",
2218 simulator_.gridView().comm());
2221 changed_well_to_group = comm.sum(
static_cast<int>(changed_well_to_group));
2222 if (changed_well_to_group) {
2223 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
2224 changed_well_group =
true;
2228 bool changed_well_individual =
false;
2232 OPM_BEGIN_PARALLEL_TRY_CATCH()
2233 for (const auto& well : well_container_) {
2234 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
2235 const bool changed_well = well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
2237 changed_well_individual = changed_well || changed_well_individual;
2240 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel: updating well controls failed: ",
2241 simulator_.gridView().comm());
2244 changed_well_individual = comm.sum(
static_cast<int>(changed_well_individual));
2245 if (changed_well_individual) {
2246 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
2247 changed_well_group =
true;
2251 const Group& fieldGroup = this->schedule().getGroup(
"FIELD", episodeIdx);
2252 this->updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState());
2254 return { changed_well_group, more_network_update };
2257 template<
typename TypeTag>
2259 BlackoilWellModel<TypeTag>::
2260 updateWellControlsDomain(DeferredLogger& deferred_logger,
const Domain& domain)
2262 if ( !this->wellsActive() ) return ;
2268 for (
const auto& well : well_container_) {
2269 if (well_domain_.at(well->name()) == domain.index) {
2270 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
2271 well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
2280 template <
typename TypeTag>
2282 BlackoilWellModel<TypeTag>::
2283 initializeWBPCalculationService()
2285 this->wbpCalcMap_.clear();
2286 this->wbpCalcMap_.resize(this->wells_ecl_.size());
2288 this->registerOpenWellsForWBPCalculation();
2290 auto wellID = std::size_t{0};
2291 for (
const auto& well : this->wells_ecl_) {
2292 this->wbpCalcMap_[wellID].wbpCalcIdx_ = this->wbpCalculationService_
2293 .createCalculator(well,
2294 this->local_parallel_well_info_[wellID],
2295 this->conn_idx_map_[wellID].local(),
2296 this->makeWellSourceEvaluatorFactory(wellID));
2301 this->wbpCalculationService_.defineCommunication();
2308 template <
typename TypeTag>
2309 data::WellBlockAveragePressures
2310 BlackoilWellModel<TypeTag>::
2311 computeWellBlockAveragePressures()
const
2313 auto wbpResult = data::WellBlockAveragePressures{};
2315 using Calculated =
typename PAvgCalculator<Scalar>::Result::WBPMode;
2316 using Output = data::WellBlockAvgPress::Quantity;
2318 this->wbpCalculationService_.collectDynamicValues();
2320 const auto numWells = this->wells_ecl_.size();
2321 for (
auto wellID = 0*numWells; wellID < numWells; ++wellID) {
2322 const auto calcIdx = this->wbpCalcMap_[wellID].wbpCalcIdx_;
2323 const auto& well = this->wells_ecl_[wellID];
2325 if (! well.hasRefDepth()) {
2331 this->wbpCalculationService_
2332 .inferBlockAveragePressures(calcIdx, well.pavg(),
2334 well.getWPaveRefDepth());
2336 const auto& result = this->wbpCalculationService_
2337 .averagePressures(calcIdx);
2339 auto& reported = wbpResult.values[well.name()];
2341 reported[Output::WBP] = result.value(Calculated::WBP);
2342 reported[Output::WBP4] = result.value(Calculated::WBP4);
2343 reported[Output::WBP5] = result.value(Calculated::WBP5);
2344 reported[Output::WBP9] = result.value(Calculated::WBP9);
2354 template <
typename TypeTag>
2355 typename ParallelWBPCalculation<typename BlackoilWellModel<TypeTag>::Scalar>::EvaluatorFactory
2356 BlackoilWellModel<TypeTag>::
2357 makeWellSourceEvaluatorFactory(
const std::vector<Well>::size_type wellIdx)
const
2359 using Span =
typename PAvgDynamicSourceData<Scalar>::template SourceDataSpan<Scalar>;
2360 using Item =
typename Span::Item;
2362 return [wellIdx,
this]() ->
typename ParallelWBPCalculation<Scalar>::Evaluator
2364 if (! this->wbpCalcMap_[wellIdx].openWellIdx_.has_value()) {
2366 return []([[maybe_unused]]
const int connIdx, Span sourceTerm)
2372 .set(Item::Pressure , 0.0)
2373 .set(Item::PoreVol , 0.0)
2374 .set(Item::MixtureDensity, 0.0);
2379 return [
this, wellPtr = this->well_container_[*this->wbpCalcMap_[wellIdx].openWellIdx_].get()]
2380 (
const int connIdx, Span sourceTerm)
2386 const auto& connIdxMap =
2387 this->conn_idx_map_[wellPtr->indexOfWell()];
2389 const auto rho = wellPtr->
2390 connectionDensity(connIdxMap.global(connIdx),
2391 connIdxMap.open(connIdx));
2394 .set(Item::Pressure , 0.0)
2395 .set(Item::PoreVol , 0.0)
2396 .set(Item::MixtureDensity, rho);
2405 template <
typename TypeTag>
2407 BlackoilWellModel<TypeTag>::
2408 registerOpenWellsForWBPCalculation()
2410 assert (this->wbpCalcMap_.size() == this->wells_ecl_.size());
2412 for (
auto& wbpCalc : this->wbpCalcMap_) {
2413 wbpCalc.openWellIdx_.reset();
2416 auto openWellIdx =
typename std::vector<WellInterfacePtr>::size_type{0};
2417 for (
const auto* openWell : this->well_container_generic_) {
2418 this->wbpCalcMap_[openWell->indexOfWell()].openWellIdx_ = openWellIdx++;
2426 template<
typename TypeTag>
2428 BlackoilWellModel<TypeTag>::
2429 updateAndCommunicate(
const int reportStepIdx,
2430 const int iterationIdx,
2431 DeferredLogger& deferred_logger)
2433 this->updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
2437 OPM_BEGIN_PARALLEL_TRY_CATCH()
2439 for (const auto& well : well_container_) {
2440 well->updateWellStateWithTarget(simulator_, this->groupState(),
2441 this->wellState(), deferred_logger);
2443 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel::updateAndCommunicate failed: ",
2444 simulator_.gridView().comm())
2445 this->updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
2448 template<typename TypeTag>
2450 BlackoilWellModel<TypeTag>::
2451 updateGroupControls(const Group& group,
2452 DeferredLogger& deferred_logger,
2453 const
int reportStepIdx,
2454 const
int iterationIdx)
2456 bool changed =
false;
2457 bool changed_hc = this->checkGroupHigherConstraints( group, deferred_logger, reportStepIdx);
2460 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
2463 bool changed_individual =
2464 BlackoilWellModelConstraints(*this).
2465 updateGroupIndividualControl(group,
2467 this->switched_inj_groups_,
2468 this->switched_prod_groups_,
2469 this->closed_offending_wells_,
2474 if (changed_individual) {
2476 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
2479 for (
const std::string& groupName : group.groups()) {
2480 bool changed_this = updateGroupControls( this->schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx,iterationIdx);
2481 changed = changed || changed_this;
2486 template<
typename TypeTag>
2492 for (
const auto& well : well_container_) {
2493 const auto& wname = well->name();
2494 const auto wasClosed = wellTestState.well_is_closed(wname);
2495 well->checkWellOperability(simulator_, this->wellState(), local_deferredLogger);
2496 const bool under_zero_target = well->wellUnderZeroGroupRateTarget(this->simulator_, this->wellState(), local_deferredLogger);
2497 well->updateWellTestState(this->wellState().well(wname), simulationTime,
true, under_zero_target, wellTestState, local_deferredLogger);
2499 if (!wasClosed && wellTestState.well_is_closed(wname)) {
2500 this->closed_this_step_.insert(wname);
2504 const Opm::Parallel::Communication comm = grid().comm();
2507 for (
const auto& [group_name, to] : this->closed_offending_wells_) {
2508 if (this->hasWell(to.second) && !this->wasDynamicallyShutThisTimeStep(to.second)) {
2509 wellTestState.close_well(to.second, WellTestConfig::Reason::GROUP, simulationTime);
2510 this->updateClosedWellsThisStep(to.second);
2511 const std::string msg =
2512 fmt::format(
"Procedure on exceeding {} limit is WELL for group {}. Well {} is {}.",
2517 global_deferredLogger.info(msg);
2521 if (this->terminal_output_) {
2527 template<
typename TypeTag>
2531 std::string& exc_msg,
2532 ExceptionType::ExcEnum& exc_type,
2535 const int np = this->numPhases();
2536 std::vector<Scalar> potentials;
2537 const auto& well = well_container_[widx];
2538 std::string cur_exc_msg;
2539 auto cur_exc_type = ExceptionType::NONE;
2541 well->computeWellPotentials(simulator_, well_state_copy, potentials, deferred_logger);
2544 OPM_PARALLEL_CATCH_CLAUSE(cur_exc_type, cur_exc_msg);
2545 if (cur_exc_type != ExceptionType::NONE) {
2546 exc_msg += fmt::format(
"\nFor well {}: {}", well->name(), cur_exc_msg);
2548 exc_type = std::max(exc_type, cur_exc_type);
2552 auto& ws = this->wellState().well(well->indexOfWell());
2553 for (
int p = 0; p < np; ++p) {
2555 ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
2561 template <
typename TypeTag>
2563 BlackoilWellModel<TypeTag>::
2564 calculateProductivityIndexValues(DeferredLogger& deferred_logger)
2566 for (
const auto& wellPtr : this->well_container_) {
2567 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2575 template <
typename TypeTag>
2577 BlackoilWellModel<TypeTag>::
2578 calculateProductivityIndexValuesShutWells(
const int reportStepIdx,
2579 DeferredLogger& deferred_logger)
2586 for (
const auto& shutWell : this->local_shut_wells_) {
2587 if (!this->wells_ecl_[shutWell].hasConnections()) {
2592 auto wellPtr = this->
template createTypedWellPointer
2593 <StandardWell<TypeTag>>(shutWell, reportStepIdx);
2595 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_, this->B_avg_,
true);
2597 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2605 template <
typename TypeTag>
2607 BlackoilWellModel<TypeTag>::
2608 calculateProductivityIndexValues(
const WellInterface<TypeTag>* wellPtr,
2609 DeferredLogger& deferred_logger)
2611 wellPtr->updateProductivityIndex(this->simulator_,
2612 this->prod_index_calc_[wellPtr->indexOfWell()],
2619 template<
typename TypeTag>
2621 BlackoilWellModel<TypeTag>::
2622 prepareTimeStep(DeferredLogger& deferred_logger)
2625 const auto episodeIdx = simulator_.episodeIndex();
2626 this->updateNetworkActiveState(episodeIdx);
2630 const bool do_prestep_network_rebalance = this->needPreStepNetworkRebalance(episodeIdx);
2632 for (
const auto& well : well_container_) {
2633 auto& events = this->wellState().well(well->indexOfWell()).events;
2634 if (events.hasEvent(WellState<Scalar>::event_mask)) {
2635 well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), deferred_logger);
2636 well->updatePrimaryVariables(simulator_, this->wellState(), deferred_logger);
2637 well->initPrimaryVariablesEvaluation();
2640 events.clearEvent(WellState<Scalar>::event_mask);
2643 if (events.hasEvent(ScheduleEvents::REQUEST_OPEN_WELL)) {
2644 events.clearEvent(ScheduleEvents::REQUEST_OPEN_WELL);
2647 if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
2649 well->solveWellEquation(simulator_, this->wellState(), this->groupState(), deferred_logger);
2650 }
catch (
const std::exception& e) {
2651 const std::string msg =
"Compute initial well solution for " + well->name() +
" initially failed. Continue with the previous rates";
2652 deferred_logger.warning(
"WELL_INITIAL_SOLVE_FAILED", msg);
2657 well->resetWellOperability();
2659 updatePrimaryVariables(deferred_logger);
2662 if (do_prestep_network_rebalance) doPreStepNetworkRebalance(deferred_logger);
2665 template<
typename TypeTag>
2667 BlackoilWellModel<TypeTag>::
2668 updateAverageFormationFactor()
2670 std::vector< Scalar > B_avg(numComponents(), Scalar() );
2671 const auto& grid = simulator_.vanguard().grid();
2672 const auto& gridView = grid.leafGridView();
2673 ElementContext elemCtx(simulator_);
2675 OPM_BEGIN_PARALLEL_TRY_CATCH();
2676 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2677 elemCtx.updatePrimaryStencil(elem);
2678 elemCtx.updatePrimaryIntensiveQuantities(0);
2680 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
2681 const auto& fs = intQuants.fluidState();
2683 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
2685 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2689 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2690 auto& B = B_avg[ compIdx ];
2692 B += 1 / fs.invB(phaseIdx).value();
2694 if constexpr (has_solvent_) {
2695 auto& B = B_avg[solventSaturationIdx];
2696 B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
2699 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel::updateAverageFormationFactor() failed: ", grid.comm())
2702 grid.comm().sum(B_avg.data(), B_avg.size());
2703 for (auto& bval : B_avg)
2705 bval /= global_num_cells_;
2714 template<
typename TypeTag>
2716 BlackoilWellModel<TypeTag>::
2717 updatePrimaryVariables(DeferredLogger& deferred_logger)
2719 for (
const auto& well : well_container_) {
2720 well->updatePrimaryVariables(simulator_, this->wellState(), deferred_logger);
2724 template<
typename TypeTag>
2726 BlackoilWellModel<TypeTag>::extractLegacyCellPvtRegionIndex_()
2728 const auto& grid = simulator_.vanguard().grid();
2729 const auto& eclProblem = simulator_.problem();
2730 const unsigned numCells = grid.size(0);
2732 this->pvt_region_idx_.resize(numCells);
2733 for (
unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
2734 this->pvt_region_idx_[cellIdx] =
2735 eclProblem.pvtRegionIndex(cellIdx);
2740 template<
typename TypeTag>
2742 BlackoilWellModel<TypeTag>::numComponents()
const
2750 int numComp = this->numPhases() < 3 ? this->numPhases() : FluidSystem::numComponents;
2751 if constexpr (has_solvent_) {
2757 template<
typename TypeTag>
2759 BlackoilWellModel<TypeTag>::extractLegacyDepth_()
2761 const auto& eclProblem = simulator_.problem();
2762 depth_.resize(local_num_cells_);
2763 for (
unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
2764 depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
2768 template<
typename TypeTag>
2769 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
2770 BlackoilWellModel<TypeTag>::
2771 getWell(
const std::string& well_name)
const
2774 auto well = std::find_if(well_container_.begin(),
2775 well_container_.end(),
2776 [&well_name](
const WellInterfacePtr& elem)->bool {
2777 return elem->name() == well_name;
2780 assert(well != well_container_.end());
2785 template<
typename TypeTag>
2787 BlackoilWellModel<TypeTag>::
2788 hasWell(
const std::string& well_name)
const
2790 return std::any_of(well_container_.begin(), well_container_.end(),
2791 [&well_name](
const WellInterfacePtr& elem) ->
bool
2793 return elem->name() == well_name;
2800 template <
typename TypeTag>
2802 BlackoilWellModel<TypeTag>::
2803 reportStepIndex()
const
2805 return std::max(this->simulator_.episodeIndex(), 0);
2812 template<
typename TypeTag>
2814 BlackoilWellModel<TypeTag>::
2815 calcResvCoeff(
const int fipnum,
2817 const std::vector<Scalar>& production_rates,
2818 std::vector<Scalar>& resv_coeff)
2820 rateConverter_->calcCoeff(fipnum, pvtreg, production_rates, resv_coeff);
2823 template<
typename TypeTag>
2825 BlackoilWellModel<TypeTag>::
2826 calcInjResvCoeff(
const int fipnum,
2828 std::vector<Scalar>& resv_coeff)
2830 rateConverter_->calcInjCoeff(fipnum, pvtreg, resv_coeff);
2834 template <
typename TypeTag>
2836 BlackoilWellModel<TypeTag>::
2837 computeWellTemperature()
2842 int np = this->numPhases();
2843 Scalar cellInternalEnergy;
2846 Scalar perfPhaseRate;
2847 const int nw = this->numLocalWells();
2848 for (
auto wellID = 0*nw; wellID < nw; ++wellID) {
2849 const Well& well = this->wells_ecl_[wellID];
2850 auto& ws = this->wellState().well(wellID);
2851 if (well.isInjector()){
2852 if( !(ws.status == WellStatus::STOP)){
2853 this->wellState().well(wellID).temperature = well.inj_temperature();
2858 std::array<Scalar,2> weighted{0.0,0.0};
2859 auto& [weighted_temperature, total_weight] = weighted;
2861 auto& well_info = this->local_parallel_well_info_[wellID].get();
2862 auto& perf_data = ws.perf_data;
2863 auto& perf_phase_rate = perf_data.phase_rates;
2865 using int_type =
decltype(this->well_perf_data_[wellID].size());
2866 for (int_type perf = 0, end_perf = this->well_perf_data_[wellID].size(); perf < end_perf; ++perf) {
2867 const int cell_idx = this->well_perf_data_[wellID][perf].cell_index;
2868 const auto& intQuants = simulator_.model().intensiveQuantities(cell_idx, 0);
2869 const auto& fs = intQuants.fluidState();
2872 Scalar cellTemperatures = fs.temperature(0).value();
2874 Scalar weight_factor = 0.0;
2875 for (
unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
2877 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2880 cellInternalEnergy = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
2881 cellBinv = fs.invB(phaseIdx).value();
2882 cellDensity = fs.density(phaseIdx).value();
2883 perfPhaseRate = perf_phase_rate[ perf*np + phaseIdx ];
2884 weight_factor += cellDensity * perfPhaseRate/cellBinv * cellInternalEnergy/cellTemperatures;
2886 weight_factor = std::abs(weight_factor)+1e-13;
2887 total_weight += weight_factor;
2888 weighted_temperature += weight_factor * cellTemperatures;
2890 well_info.communication().sum(weighted.data(), 2);
2891 this->wellState().well(wellID).temperature = weighted_temperature/total_weight;
2896 template <
typename TypeTag>
2898 BlackoilWellModel<TypeTag>::
2899 logPrimaryVars()
const
2901 std::ostringstream os;
2902 for (
const auto& w : well_container_) {
2903 os << w->name() <<
":";
2904 auto pv = w->getPrimaryVars();
2905 for (
const Scalar v : pv) {
2910 OpmLog::debug(os.str());
2915 template <
typename TypeTag>
2916 std::vector<typename BlackoilWellModel<TypeTag>::Scalar>
2917 BlackoilWellModel<TypeTag>::
2918 getPrimaryVarsDomain(
const Domain& domain)
const
2920 std::vector<Scalar> ret;
2921 for (
const auto& well : well_container_) {
2922 if (well_domain_.at(well->name()) == domain.index) {
2923 const auto& pv = well->getPrimaryVars();
2924 ret.insert(ret.end(), pv.begin(), pv.end());
2932 template <
typename TypeTag>
2934 BlackoilWellModel<TypeTag>::
2935 setPrimaryVarsDomain(
const Domain& domain,
const std::vector<Scalar>& vars)
2937 std::size_t offset = 0;
2938 for (
auto& well : well_container_) {
2939 if (well_domain_.at(well->name()) == domain.index) {
2940 int num_pri_vars = well->setPrimaryVars(vars.begin() + offset);
2941 offset += num_pri_vars;
2944 assert(offset == vars.size());
2949 template <
typename TypeTag>
2951 BlackoilWellModel<TypeTag>::
2952 setupDomains(
const std::vector<Domain>& domains)
2954 OPM_BEGIN_PARALLEL_TRY_CATCH();
2959 for (
const auto& wellPtr : this->well_container_) {
2960 const int first_well_cell = wellPtr->cells().front();
2961 for (
const auto& domain : domains) {
2962 auto cell_present = [&domain](
const auto cell)
2964 return std::binary_search(domain.cells.begin(),
2965 domain.cells.end(), cell);
2968 if (cell_present(first_well_cell)) {
2971 well_domain_[wellPtr->name()] = domain.index;
2974 for (
int well_cell : wellPtr->cells()) {
2975 if (! cell_present(well_cell)) {
2976 OPM_THROW(std::runtime_error,
2977 fmt::format(
"Well '{}' found on multiple domains.",
2984 OPM_END_PARALLEL_TRY_CATCH(
"BlackoilWellModel::setupDomains(): well found on multiple domains.",
2985 simulator_.gridView().comm());
2988 const Opm::Parallel::Communication& comm = grid().comm();
2989 const int rank = comm.rank();
2990 DeferredLogger local_log;
2991 if (!well_domain_.empty()) {
2992 std::ostringstream os;
2993 os <<
"Well name Rank Domain\n";
2994 for (
const auto& [wname, domain] : well_domain_) {
2995 os << wname << std::setw(19 - wname.size()) << rank << std::setw(12) << domain <<
'\n';
2997 local_log.debug(os.str());
3000 if (this->terminal_output_) {
3001 global_log.logMessages();
3005 well_local_cells_.clear();
3006 well_local_cells_.reserve(well_container_.size(), 10);
3007 std::vector<int> local_cells;
3008 for (
const auto& well : well_container_) {
3009 const auto& global_cells = well->cells();
3010 const int domain_index = well_domain_.at(well->name());
3011 const auto& domain_cells = domains[domain_index].cells;
3012 local_cells.resize(global_cells.size());
3016 for (
size_t i = 0; i < global_cells.size(); ++i) {
3017 auto it = std::lower_bound(domain_cells.begin(), domain_cells.end(), global_cells[i]);
3018 if (it != domain_cells.end() && *it == global_cells[i]) {
3019 local_cells[i] = std::distance(domain_cells.begin(), it);
3021 OPM_THROW(std::runtime_error, fmt::format(
"Cell {} not found in domain {}",
3022 global_cells[i], domain_index));
3025 well_local_cells_.appendRow(local_cells.begin(), local_cells.end());
Class for handling the blackoil well model.
Definition WellInterface.hpp:31
Definition DeferredLogger.hpp:57
void logMessages()
Log all messages to the OpmLog backends, and clear the message container.
Definition DeferredLogger.cpp:85
Class for serializing and broadcasting data using MPI.
Definition MPISerializer.hpp:31
void broadcast(T &data, int root=0)
Serialize and broadcast on root process, de-serialize on others.
Definition MPISerializer.hpp:46
Definition WellGroupHelpers.hpp:48
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:62
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
ConvergenceReport gatherConvergenceReport(const ConvergenceReport &local_report, Parallel::Communication mpi_communicator)
Create a global convergence report combining local (per-process) reports.
Definition gatherConvergenceReport.cpp:171
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition phaseUsageFromDeck.cpp:137
Opm::DeferredLogger gatherDeferredLogger(const Opm::DeferredLogger &local_deferredlogger, Opm::Parallel::Communication)
Create a global log combining local logs.
Definition gatherDeferredLogger.cpp:168