My Project
BlackoilWellModel_impl.hpp
1 /*
2  Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics.
3  Copyright 2016 - 2018 Equinor ASA.
4  Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
5  Copyright 2016 - 2018 Norce AS
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24 #include <opm/core/props/phaseUsageFromDeck.hpp>
25 #include <opm/grid/utility/cartesianToCompressed.hpp>
26 #include <opm/input/eclipse/Units/UnitSystem.hpp>
27 
28 #include <opm/simulators/wells/VFPProperties.hpp>
29 #include <opm/simulators/utils/MPIPacker.hpp>
30 #include <opm/simulators/linalg/bda/WellContributions.hpp>
31 
32 #if HAVE_MPI
33 #include <ebos/eclmpiserializer.hh>
34 #endif
35 
36 #include <algorithm>
37 #include <utility>
38 
39 #include <fmt/format.h>
40 
41 namespace Opm {
42  template<typename TypeTag>
43  BlackoilWellModel<TypeTag>::
44  BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage)
45  : BlackoilWellModelGeneric(ebosSimulator.vanguard().schedule(),
46  ebosSimulator.vanguard().summaryState(),
47  ebosSimulator.vanguard().eclState(),
48  phase_usage,
49  ebosSimulator.gridView().comm())
50  , ebosSimulator_(ebosSimulator)
51  {
52  terminal_output_ = ((ebosSimulator.gridView().comm().rank() == 0) &&
53  EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput));
54 
55  local_num_cells_ = ebosSimulator_.gridView().size(0);
56 
57  // Number of cells the global grid view
58  global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();
59 
60  {
61  auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
62 
63  this->parallel_well_info_.reserve(parallel_wells.size());
64  for( const auto& name_bool : parallel_wells) {
65  this->parallel_well_info_.emplace_back(name_bool, grid().comm());
66  }
67  }
68 
69  this->alternative_well_rate_init_ =
70  EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
71  }
72 
73  template<typename TypeTag>
74  BlackoilWellModel<TypeTag>::
75  BlackoilWellModel(Simulator& ebosSimulator) :
76  BlackoilWellModel(ebosSimulator, phaseUsageFromDeck(ebosSimulator.vanguard().eclState()))
77  {}
78 
79 
80  template<typename TypeTag>
81  void
82  BlackoilWellModel<TypeTag>::
83  init()
84  {
85  extractLegacyCellPvtRegionIndex_();
86  extractLegacyDepth_();
87 
88  gravity_ = ebosSimulator_.problem().gravity()[2];
89 
90  initial_step_ = true;
91 
92  // add the eWoms auxiliary module for the wells to the list
93  ebosSimulator_.model().addAuxiliaryModule(this);
94 
95  is_cell_perforated_.resize(local_num_cells_, false);
96  }
97 
98 
99  template<typename TypeTag>
100  void
101  BlackoilWellModel<TypeTag>::
102  initWellContainer(const int reportStepIdx)
103  {
104  const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
105  + ScheduleEvents::NEW_WELL;
106  const auto& events = schedule()[reportStepIdx].wellgroup_events();
107  for (auto& wellPtr : this->well_container_) {
108  const bool well_opened_this_step = report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask);
109  wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
110  this->local_num_cells_, this->B_avg_, well_opened_this_step);
111  }
112  }
113 
114  template<typename TypeTag>
115  void
116  BlackoilWellModel<TypeTag>::
117  addNeighbors(std::vector<NeighborSet>& neighbors) const
118  {
119  if (!param_.matrix_add_well_contributions_) {
120  return;
121  }
122 
123  // Create cartesian to compressed mapping
124  const auto& schedule_wells = schedule().getWellsatEnd();
125 
126  // initialize the additional cell connections introduced by wells.
127  for (const auto& well : schedule_wells)
128  {
129  std::vector<int> wellCells;
130  // All possible connections of the well
131  const auto& connectionSet = well.getConnections();
132  wellCells.reserve(connectionSet.size());
133 
134  for ( size_t c=0; c < connectionSet.size(); c++ )
135  {
136  const auto& connection = connectionSet.get(c);
137  int compressed_idx = compressedIndexForInterior(connection.global_index());
138  if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
139  wellCells.push_back(compressed_idx);
140  }
141  }
142 
143  for (int cellIdx : wellCells) {
144  neighbors[cellIdx].insert(wellCells.begin(),
145  wellCells.end());
146  }
147  }
148  }
149 
150  template<typename TypeTag>
151  void
152  BlackoilWellModel<TypeTag>::
153  linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
154  {
155  if (!param_.matrix_add_well_contributions_)
156  {
157  OPM_BEGIN_PARALLEL_TRY_CATCH();
158  {
159  // if the well contributions are not supposed to be included explicitly in
160  // the matrix, we only apply the vector part of the Schur complement here.
161  for (const auto& well: well_container_) {
162  // r = r - duneC_^T * invDuneD_ * resWell_
163  well->apply(res);
164  }
165  }
166  OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ",
167  ebosSimulator_.gridView().comm());
168  return;
169  }
170 
171  for (const auto& well: well_container_) {
172  well->addWellContributions(jacobian);
173 
174  // applying the well residual to reservoir residuals
175  // r = r - duneC_^T * invDuneD_ * resWell_
176  well->apply(res);
177  }
178  }
179 
180 
181  template<typename TypeTag>
182  void
183  BlackoilWellModel<TypeTag>::
184  beginReportStep(const int timeStepIdx)
185  {
186  DeferredLogger local_deferredLogger;
187 
188  report_step_starts_ = true;
189 
190  const Grid& grid = ebosSimulator_.vanguard().grid();
191  const auto& summaryState = ebosSimulator_.vanguard().summaryState();
192  // Make wells_ecl_ contain only this partition's wells.
193  wells_ecl_ = getLocalWells(timeStepIdx);
194  this->local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
195 
196  // at least initializeWellState might be throw
197  // exception in opm-material (UniformTabulated2DFunction.hpp)
198  // playing it safe by extending the scope a bit.
199  OPM_BEGIN_PARALLEL_TRY_CATCH();
200  {
201 
202  // The well state initialize bhp with the cell pressure in the top cell.
203  // We must therefore provide it with updated cell pressures
204  this->initializeWellPerfData();
205  this->initializeWellState(timeStepIdx, summaryState);
206 
207  // handling MS well related
208  if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
209  this->wellState().initWellStateMSWell(wells_ecl_, &this->prevWellState());
210  }
211 
212  const Group& fieldGroup = schedule().getGroup("FIELD", timeStepIdx);
213  WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState());
214 
215  // Compute reservoir volumes for RESV controls.
216  rateConverter_.reset(new RateConverterType (phase_usage_,
217  std::vector<int>(local_num_cells_, 0)));
218  rateConverter_->template defineState<ElementContext>(ebosSimulator_);
219 
220  // Compute regional average pressures used by gpmaint
221  if (schedule_[timeStepIdx].has_gpmaint()) {
222  const auto& fp = this->eclState_.fieldProps();
223  const auto& fipnum = fp.get_int("FIPNUM");
224  regionalAveragePressureCalculator_.reset(new AverageRegionalPressureType (phase_usage_,fipnum));
225  }
226 
227  {
228  const auto& sched_state = this->schedule()[timeStepIdx];
229  // update VFP properties
230  vfp_properties_.reset(new VFPProperties( sched_state.vfpinj(), sched_state.vfpprod(), this->prevWellState()));
231  this->initializeWellProdIndCalculators();
232  if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) {
233  this->runWellPIScaling(timeStepIdx, local_deferredLogger);
234  }
235  }
236  }
237  OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "beginReportStep() failed: ",
238  terminal_output_, grid.comm());
239  // Store the current well state, to be able to recover in the case of failed iterations
240  this->commitWGState();
241  }
242 
243 
244  // called at the beginning of a time step
245  template<typename TypeTag>
246  void
247  BlackoilWellModel<TypeTag>::
248  beginTimeStep()
249  {
250  updatePerforationIntensiveQuantities();
251  updateAverageFormationFactor();
252  DeferredLogger local_deferredLogger;
253  switched_prod_groups_.clear();
254  switched_inj_groups_.clear();
255 
256  this->resetWGState();
257  const int reportStepIdx = ebosSimulator_.episodeIndex();
258  updateAndCommunicateGroupData(reportStepIdx,
259  ebosSimulator_.model().newtonMethod().numIterations());
260  this->wellState().gliftTimeStepInit();
261  const double simulationTime = ebosSimulator_.time();
262  OPM_BEGIN_PARALLEL_TRY_CATCH();
263  {
264  // test wells
265  wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
266 
267  // create the well container
268  createWellContainer(reportStepIdx);
269 
270  // Wells are active if they are active wells on at least one process.
271  const Grid& grid = ebosSimulator_.vanguard().grid();
272  wells_active_ = !this->well_container_.empty();
273  wells_active_ = grid.comm().max(wells_active_);
274 
275  // do the initialization for all the wells
276  // TODO: to see whether we can postpone of the intialization of the well containers to
277  // optimize the usage of the following several member variables
278  this->initWellContainer(reportStepIdx);
279 
280  // update the updated cell flag
281  std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
282  for (auto& well : well_container_) {
283  well->updatePerforatedCell(is_cell_perforated_);
284  }
285 
286  // calculate the efficiency factors for each well
287  calculateEfficiencyFactors(reportStepIdx);
288 
289  if constexpr (has_polymer_)
290  {
291  if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
292  setRepRadiusPerfLength();
293  }
294  }
295 
296  }
297  OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "beginTimeStep() failed: ",
298  terminal_output_, ebosSimulator_.vanguard().grid().comm());
299 
300  for (auto& well : well_container_) {
301  well->setVFPProperties(vfp_properties_.get());
302  well->setGuideRate(&guideRate_);
303  }
304 
305  // Close completions due to economical reasons
306  for (auto& well : well_container_) {
307  well->closeCompletions(wellTestState());
308  }
309 
310  if (alternative_well_rate_init_) {
311  // Update the well rates of well_state_, if only single-phase rates, to
312  // have proper multi-phase rates proportional to rates at bhp zero.
313  // This is done only for producers, as injectors will only have a single
314  // nonzero phase anyway.
315  for (auto& well : well_container_) {
316  if (well->isProducer()) {
317  well->updateWellStateRates(ebosSimulator_, this->wellState(), local_deferredLogger);
318  }
319  }
320  }
321 
322  // calculate the well potentials
323  try {
324  updateWellPotentials(reportStepIdx,
325  /*onlyAfterEvent*/true,
326  ebosSimulator_.vanguard().summaryConfig(),
327  local_deferredLogger);
328  } catch ( std::runtime_error& e ) {
329  const std::string msg = "A zero well potential is returned for output purposes. ";
330  local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
331  }
332 
333  //update guide rates
334  const auto& comm = ebosSimulator_.vanguard().grid().comm();
335  const auto& summaryState = ebosSimulator_.vanguard().summaryState();
336  std::vector<double> pot(numPhases(), 0.0);
337  const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
338  WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
339  this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
340  std::string exc_msg;
341  auto exc_type = ExceptionType::NONE;
342  // update gpmaint targets
343  if (schedule_[reportStepIdx].has_gpmaint()) {
344  regionalAveragePressureCalculator_->template defineState<ElementContext>(ebosSimulator_);
345  const double dt = ebosSimulator_.timeStepSize();
346  WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup,
347  schedule_, *regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState());
348  }
349  try {
350  // Compute initial well solution for new wells and injectors that change injection type i.e. WAG.
351  for (auto& well : well_container_) {
352  const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
353  + ScheduleEvents::INJECTION_TYPE_CHANGED
354  + ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
355  + ScheduleEvents::NEW_WELL;
356 
357  const auto& events = schedule()[reportStepIdx].wellgroup_events();
358  const bool event = report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
359  const bool dyn_status_change = this->wellState().well(well->name()).status
360  != this->prevWellState().well(well->name()).status;
361 
362  if (event || dyn_status_change) {
363  try {
364  well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), local_deferredLogger);
365  well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), local_deferredLogger);
366  well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), local_deferredLogger);
367  } catch (const std::exception& e) {
368  const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates";
369  local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
370  }
371  }
372  }
373  }
374  // Catch clauses for all errors setting exc_type and exc_msg
375  OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
376 
377  if (exc_type != ExceptionType::NONE) {
378  const std::string msg = "Compute initial well solution for new wells failed. Continue with zero initial rates";
379  local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
380  }
381 
382  logAndCheckForExceptionsAndThrow(local_deferredLogger,
383  exc_type, "beginTimeStep() failed: " + exc_msg, terminal_output_, comm);
384 
385  }
386 
387  template<typename TypeTag>
388  void
389  BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx,
390  const double simulationTime,
391  DeferredLogger& deferred_logger)
392  {
393  const auto& wtest_config = schedule()[timeStepIdx].wtest_config();
394  if (!wtest_config.empty()) { // there is a WTEST request
395  const std::vector<std::string> wellsForTesting = wellTestState()
396  .test_wells(wtest_config, simulationTime);
397  for (const std::string& well_name : wellsForTesting) {
398 
399  const Well& wellEcl = schedule().getWell(well_name, timeStepIdx);
400  if (wellEcl.getStatus() == Well::Status::SHUT)
401  continue;
402 
403  WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
404  // some preparation before the well can be used
405  well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_, true);
406 
407  double well_efficiency_factor = wellEcl.getEfficiencyFactor();
408  WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx),
409  schedule(), timeStepIdx, well_efficiency_factor);
410 
411  well->setWellEfficiencyFactor(well_efficiency_factor);
412  well->setVFPProperties(vfp_properties_.get());
413  well->setGuideRate(&guideRate_);
414 
415  well->wellTesting(ebosSimulator_, simulationTime, this->wellState(), this->groupState(), wellTestState(), deferred_logger);
416  }
417  }
418  }
419 
420 
421 
422 
423 
424  // called at the end of a report step
425  template<typename TypeTag>
426  void
427  BlackoilWellModel<TypeTag>::
428  endReportStep()
429  {
430  // Clear the communication data structures for above values.
431  for (auto&& pinfo : this->local_parallel_well_info_)
432  {
433  pinfo.get().clear();
434  }
435  }
436 
437 
438 
439 
440 
441  // called at the end of a report step
442  template<typename TypeTag>
443  const SimulatorReportSingle&
444  BlackoilWellModel<TypeTag>::
445  lastReport() const {return last_report_; }
446 
447 
448 
449 
450 
451  // called at the end of a time step
452  template<typename TypeTag>
453  void
454  BlackoilWellModel<TypeTag>::
455  timeStepSucceeded(const double& simulationTime, const double dt)
456  {
457  this->closed_this_step_.clear();
458 
459  // time step is finished and we are not any more at the beginning of an report step
460  report_step_starts_ = false;
461  const int reportStepIdx = ebosSimulator_.episodeIndex();
462 
463  DeferredLogger local_deferredLogger;
464  for (const auto& well : well_container_) {
465  if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
466  well->updateWaterThroughput(dt, this->wellState());
467  }
468  }
469  // report well switching
470  for (const auto& well : well_container_) {
471  well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);
472  }
473  // report group switching
474  if (terminal_output_) {
475 
476  for (const auto& [name, to] : switched_prod_groups_) {
477  const Group::ProductionCMode& oldControl = this->prevWGState().group_state.production_control(name);
478  std::string from = Group::ProductionCMode2String(oldControl);
479  if (to != from) {
480  std::string msg = " Production Group " + name
481  + " control mode changed from ";
482  msg += from;
483  msg += " to " + to;
484  local_deferredLogger.info(msg);
485  }
486  }
487  for (const auto& [key, to] : switched_inj_groups_) {
488  const std::string& name = key.first;
489  const Opm::Phase& phase = key.second;
490 
491  const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(name, phase);
492  std::string from = Group::InjectionCMode2String(oldControl);
493  if (to != from) {
494  std::string msg = " Injection Group " + name
495  + " control mode changed from ";
496  msg += from;
497  msg += " to " + to;
498  local_deferredLogger.info(msg);
499  }
500  }
501  }
502 
503  // update the rate converter with current averages pressures etc in
504  rateConverter_->template defineState<ElementContext>(ebosSimulator_);
505 
506  // calculate the well potentials
507  try {
508  updateWellPotentials(reportStepIdx,
509  /*onlyAfterEvent*/false,
510  ebosSimulator_.vanguard().summaryConfig(),
511  local_deferredLogger);
512  } catch ( std::runtime_error& e ) {
513  const std::string msg = "A zero well potential is returned for output purposes. ";
514  local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
515  }
516 
517  updateWellTestState(simulationTime, wellTestState());
518 
519  // check group sales limits at the end of the timestep
520  const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx);
521  checkGconsaleLimits(fieldGroup, this->wellState(),
522  ebosSimulator_.episodeIndex(), local_deferredLogger);
523 
524  this->calculateProductivityIndexValues(local_deferredLogger);
525 
526  this->commitWGState();
527 
528  const Opm::Parallel::Communication& comm = grid().comm();
529  DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
530  if (terminal_output_) {
531  global_deferredLogger.logMessages();
532  }
533 
534  //reporting output temperatures
535  this->computeWellTemperature();
536  }
537 
538 
539  template<typename TypeTag>
540  void
541  BlackoilWellModel<TypeTag>::
542  computeTotalRatesForDof(RateVector& rate,
543  unsigned elemIdx) const
544  {
545  rate = 0;
546 
547  if (!is_cell_perforated_[elemIdx])
548  return;
549 
550  for (const auto& well : well_container_)
551  well->addCellRates(rate, elemIdx);
552  }
553 
554 
555  template<typename TypeTag>
556  template <class Context>
557  void
558  BlackoilWellModel<TypeTag>::
559  computeTotalRatesForDof(RateVector& rate,
560  const Context& context,
561  unsigned spaceIdx,
562  unsigned timeIdx) const
563  {
564  rate = 0;
565  int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
566 
567  if (!is_cell_perforated_[elemIdx])
568  return;
569 
570  for (const auto& well : well_container_)
571  well->addCellRates(rate, elemIdx);
572  }
573 
574 
575 
576  template<typename TypeTag>
577  void
578  BlackoilWellModel<TypeTag>::
579  initializeWellState(const int timeStepIdx,
580  const SummaryState& summaryState)
581  {
582  std::vector<double> cellPressures(this->local_num_cells_, 0.0);
583  ElementContext elemCtx(ebosSimulator_);
584 
585  const auto& gridView = ebosSimulator_.vanguard().gridView();
586 
587  OPM_BEGIN_PARALLEL_TRY_CATCH();
588  for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
589  elemCtx.updatePrimaryStencil(elem);
590  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
591 
592  const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState();
593  // copy of get perfpressure in Standard well except for value
594  double& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)];
595  if (Indices::oilEnabled) {
596  perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
597  } else if (Indices::waterEnabled) {
598  perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
599  } else {
600  perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
601  }
602  }
603  OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ", ebosSimulator_.vanguard().grid().comm());
604 
605  this->wellState().init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx,
606  &this->prevWellState(), well_perf_data_,
607  summaryState);
608  }
609 
610 
611 
612 
613 
614  template<typename TypeTag>
615  void
616  BlackoilWellModel<TypeTag>::
617  createWellContainer(const int time_step)
618  {
619  DeferredLogger local_deferredLogger;
620 
621  const int nw = numLocalWells();
622 
623  well_container_.clear();
624 
625  if (nw > 0) {
626  well_container_.reserve(nw);
627 
628  for (int w = 0; w < nw; ++w) {
629  const Well& well_ecl = wells_ecl_[w];
630 
631  if (well_ecl.getConnections().empty()) {
632  // No connections in this well. Nothing to do.
633  continue;
634  }
635 
636  const std::string& well_name = well_ecl.name();
637  const auto well_status = this->schedule()
638  .getWell(well_name, time_step).getStatus();
639 
640  if ((well_ecl.getStatus() == Well::Status::SHUT) ||
641  (well_status == Well::Status::SHUT))
642  {
643  // Due to ACTIONX the well might have been closed behind our back.
644  if (well_ecl.getStatus() != Well::Status::SHUT) {
645  this->closed_this_step_.insert(well_name);
646  this->wellState().shutWell(w);
647  }
648 
649  continue;
650  }
651 
652  // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
653  if (this->wellTestState().well_is_closed(well_name)) {
654  // TODO: more checking here, to make sure this standard more specific and complete
655  // maybe there is some WCON keywords will not open the well
656  auto& events = this->wellState().well(w).events;
657  if (events.hasEvent(WellState::event_mask)) {
658  if (wellTestState().lastTestTime(well_name) == ebosSimulator_.time()) {
659  // The well was shut this timestep, we are most likely retrying
660  // a timestep without the well in question, after it caused
661  // repeated timestep cuts. It should therefore not be opened,
662  // even if it was new or received new targets this report step.
663  events.clearEvent(WellState::event_mask);
664  } else {
665  wellTestState().open_well(well_name);
666  wellTestState().open_completions(well_name);
667  }
668  }
669  }
670 
671  // TODO: should we do this for all kinds of closing reasons?
672  // something like wellTestState().hasWell(well_name)?
673  bool wellIsStopped = false;
674  if (wellTestState().well_is_closed(well_name))
675  {
676  if (well_ecl.getAutomaticShutIn()) {
677  // shut wells are not added to the well container
678  this->wellState().shutWell(w);
679  continue;
680  } else {
681  if (!well_ecl.getAllowCrossFlow()) {
682  // stopped wells where cross flow is not allowed
683  // are not added to the well container
684  this->wellState().shutWell(w);
685  continue;
686  }
687  // stopped wells are added to the container but marked as stopped
688  this->wellState().stopWell(w);
689  wellIsStopped = true;
690  }
691  }
692 
693  // If a production well disallows crossflow and its
694  // (prediction type) rate control is zero, then it is effectively shut.
695  if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) {
696  const auto& summaryState = ebosSimulator_.vanguard().summaryState();
697  const auto prod_controls = well_ecl.productionControls(summaryState);
698 
699  auto is_zero = [](const double x)
700  {
701  return std::isfinite(x) && !std::isnormal(x);
702  };
703 
704  bool zero_rate_control = false;
705  switch (prod_controls.cmode) {
706  case Well::ProducerCMode::ORAT:
707  zero_rate_control = is_zero(prod_controls.oil_rate);
708  break;
709 
710  case Well::ProducerCMode::WRAT:
711  zero_rate_control = is_zero(prod_controls.water_rate);
712  break;
713 
714  case Well::ProducerCMode::GRAT:
715  zero_rate_control = is_zero(prod_controls.gas_rate);
716  break;
717 
718  case Well::ProducerCMode::LRAT:
719  zero_rate_control = is_zero(prod_controls.liquid_rate);
720  break;
721 
722  case Well::ProducerCMode::RESV:
723  zero_rate_control = is_zero(prod_controls.resv_rate);
724  break;
725  default:
726  // Might still have zero rate controls, but is pressure controlled.
727  zero_rate_control = false;
728  break;
729  }
730 
731  if (zero_rate_control) {
732  // Treat as shut, do not add to container.
733  local_deferredLogger.info(" Well shut due to zero rate control and disallowing crossflow: " + well_ecl.name());
734  this->wellState().shutWell(w);
735  continue;
736  }
737  }
738 
739  if (well_status == Well::Status::STOP) {
740  this->wellState().stopWell(w);
741  wellIsStopped = true;
742  }
743 
744  well_container_.emplace_back(this->createWellPointer(w, time_step));
745 
746  if (wellIsStopped)
747  well_container_.back()->stopWell();
748  }
749  }
750 
751  // Collect log messages and print.
752 
753  const Opm::Parallel::Communication& comm = grid().comm();
754  DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
755  if (terminal_output_) {
756  global_deferredLogger.logMessages();
757  }
758 
759  well_container_generic_.clear();
760  for (auto& w : well_container_)
761  well_container_generic_.push_back(w.get());
762  }
763 
764 
765 
766 
767 
768  template <typename TypeTag>
769  typename BlackoilWellModel<TypeTag>::WellInterfacePtr
770  BlackoilWellModel<TypeTag>::
771  createWellPointer(const int wellID, const int time_step) const
772  {
773  const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
774 
775  if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
776  return this->template createTypedWellPointer<StandardWell<TypeTag>>(wellID, time_step);
777  }
778  else {
779  return this->template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, time_step);
780  }
781  }
782 
783 
784 
785 
786 
787  template <typename TypeTag>
788  template <typename WellType>
789  std::unique_ptr<WellType>
790  BlackoilWellModel<TypeTag>::
791  createTypedWellPointer(const int wellID, const int time_step) const
792  {
793  // Use the pvtRegionIdx from the top cell
794  const auto& perf_data = this->well_perf_data_[wellID];
795 
796  // Cater for case where local part might have no perforations.
797  const auto pvtreg = perf_data.empty()
798  ? 0 : pvt_region_idx_[perf_data.front().cell_index];
799 
800  const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get();
801  const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
802 
803  return std::make_unique<WellType>(this->wells_ecl_[wellID],
804  parallel_well_info,
805  time_step,
806  this->param_,
807  *this->rateConverter_,
808  global_pvtreg,
809  this->numComponents(),
810  this->numPhases(),
811  wellID,
812  perf_data);
813  }
814 
815 
816 
817 
818 
819  template<typename TypeTag>
820  typename BlackoilWellModel<TypeTag>::WellInterfacePtr
821  BlackoilWellModel<TypeTag>::
822  createWellForWellTest(const std::string& well_name,
823  const int report_step,
824  DeferredLogger& deferred_logger) const
825  {
826  // Finding the location of the well in wells_ecl
827  const int nw_wells_ecl = wells_ecl_.size();
828  int index_well_ecl = 0;
829  for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
830  if (well_name == wells_ecl_[index_well_ecl].name()) {
831  break;
832  }
833  }
834  // It should be able to find in wells_ecl.
835  if (index_well_ecl == nw_wells_ecl) {
836  OPM_DEFLOG_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ", deferred_logger);
837  }
838 
839  return this->createWellPointer(index_well_ecl, report_step);
840  }
841 
842 
843 
844 
845 
846  template<typename TypeTag>
847  void
848  BlackoilWellModel<TypeTag>::
849  assemble(const int iterationIdx,
850  const double dt)
851  {
852 
853  DeferredLogger local_deferredLogger;
854  if (this->glift_debug) {
855  const std::string msg = fmt::format(
856  "assemble() : iteration {}" , iterationIdx);
857  gliftDebug(msg, local_deferredLogger);
858  }
859  last_report_ = SimulatorReportSingle();
860  Dune::Timer perfTimer;
861  perfTimer.start();
862 
863  if ( ! wellsActive() ) {
864  return;
865  }
866 
867  updatePerforationIntensiveQuantities();
868 
869  if (iterationIdx == 0) {
870  // try-catch is needed here as updateWellControls
871  // contains global communication and has either to
872  // be reached by all processes or all need to abort
873  // before.
874  OPM_BEGIN_PARALLEL_TRY_CATCH();
875  {
876  calculateExplicitQuantities(local_deferredLogger);
877  prepareTimeStep(local_deferredLogger);
878  }
879  OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed (It=0): ",
880  terminal_output_, grid().comm());
881  }
882 
883  const bool well_group_control_changed = assembleImpl(iterationIdx, dt, 0, local_deferredLogger);
884 
885  // if group or well control changes we don't consider the
886  // case converged
887  last_report_.well_group_control_changed = well_group_control_changed;
888  last_report_.assemble_time_well += perfTimer.stop();
889  }
890 
891 
892 
893 
894 
895  template<typename TypeTag>
896  bool
897  BlackoilWellModel<TypeTag>::
898  assembleImpl(const int iterationIdx,
899  const double dt,
900  const std::size_t recursion_level,
901  DeferredLogger& local_deferredLogger)
902  {
903 
904  auto [well_group_control_changed, network_changed, network_imbalance] = updateWellControls(local_deferredLogger);
905 
906  bool alq_updated = false;
907  OPM_BEGIN_PARALLEL_TRY_CATCH();
908  {
909  // Set the well primary variables based on the value of well solutions
910  initPrimaryVariablesEvaluation();
911 
912  alq_updated = maybeDoGasLiftOptimize(local_deferredLogger);
913  assembleWellEq(dt, local_deferredLogger);
914  }
915  OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed: ",
916  terminal_output_, grid().comm());
917 
918  //update guide rates
919  const int reportStepIdx = ebosSimulator_.episodeIndex();
920  if (alq_updated || guideRateUpdateIsNeeded(reportStepIdx)) {
921  const double simulationTime = ebosSimulator_.time();
922  const auto& comm = ebosSimulator_.vanguard().grid().comm();
923  const auto& summaryState = ebosSimulator_.vanguard().summaryState();
924  std::vector<double> pot(numPhases(), 0.0);
925  const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
926  WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
927  this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
928  }
929 
930  // Maybe do a recursive call to iterate network and well controls.
931  if (network_changed) {
932  if (shouldBalanceNetwork(reportStepIdx, iterationIdx)) {
933  const auto& balance = schedule()[reportStepIdx].network_balance();
934  // Iterate if not converged, and number of iterations is not yet max (NETBALAN item 3).
935  if (recursion_level < balance.pressure_max_iter() && network_imbalance > balance.pressure_tolerance()) {
936  well_group_control_changed = assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
937  }
938  }
939  }
940  return well_group_control_changed;
941  }
942 
943 
944 
945 
946 
947  template<typename TypeTag>
948  bool
949  BlackoilWellModel<TypeTag>::
950  maybeDoGasLiftOptimize(DeferredLogger& deferred_logger)
951  {
952  bool do_glift_optimization = false;
953  int num_wells_changed = 0;
954  const double simulation_time = ebosSimulator_.time();
955  const double min_wait = ebosSimulator_.vanguard().schedule().glo(ebosSimulator_.episodeIndex()).min_wait();
956  // We only optimize if a min_wait time has past.
957  // If all_newton is true we still want to optimize several times pr timestep
958  // i.e. we also optimize if check simulation_time == last_glift_opt_time_
959  // that is when the last_glift_opt_time is already updated with the current time step
960  if ( simulation_time == last_glift_opt_time_ || simulation_time >= (last_glift_opt_time_ + min_wait)) {
961  do_glift_optimization = true;
962  last_glift_opt_time_ = simulation_time;
963  }
964 
965  if (do_glift_optimization) {
966  GLiftOptWells glift_wells;
967  GLiftProdWells prod_wells;
968  GLiftWellStateMap state_map;
969  // NOTE: To make GasLiftGroupInfo (see below) independent of the TypeTag
970  // associated with *this (i.e. BlackoilWellModel<TypeTag>) we observe
971  // that GasLiftGroupInfo's only dependence on *this is that it needs to
972  // access the eclipse Wells in the well container (the eclipse Wells
973  // themselves are independent of the TypeTag).
974  // Hence, we extract them from the well container such that we can pass
975  // them to the GasLiftGroupInfo constructor.
976  GLiftEclWells ecl_well_map;
977  initGliftEclWellMap(ecl_well_map);
978  GasLiftGroupInfo group_info {
979  ecl_well_map,
980  ebosSimulator_.vanguard().schedule(),
981  ebosSimulator_.vanguard().summaryState(),
982  ebosSimulator_.episodeIndex(),
983  ebosSimulator_.model().newtonMethod().numIterations(),
984  phase_usage_,
985  deferred_logger,
986  this->wellState(),
987  ebosSimulator_.vanguard().grid().comm(),
988  this->glift_debug
989  };
990  group_info.initialize();
991  gasLiftOptimizationStage1(
992  deferred_logger, prod_wells, glift_wells, group_info, state_map);
993  gasLiftOptimizationStage2(
994  deferred_logger, prod_wells, glift_wells, state_map,
995  ebosSimulator_.episodeIndex());
996  if (this->glift_debug) gliftDebugShowALQ(deferred_logger);
997  num_wells_changed = glift_wells.size();
998  }
999  num_wells_changed = this->comm_.sum(num_wells_changed);
1000  return num_wells_changed > 0;
1001  }
1002 
1003  template<typename TypeTag>
1004  void
1005  BlackoilWellModel<TypeTag>::
1006  gasLiftOptimizationStage1(DeferredLogger& deferred_logger,
1007  GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
1008  GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map)
1009  {
1010  auto comm = ebosSimulator_.vanguard().grid().comm();
1011  int num_procs = comm.size();
1012  // NOTE: Gas lift optimization stage 1 seems to be difficult
1013  // to do in parallel since the wells are optimized on different
1014  // processes and each process needs to know the current ALQ allocated
1015  // to each group it is a memeber of in order to check group limits and avoid
1016  // allocating more ALQ than necessary. (Surplus ALQ is removed in
1017  // stage 2). In stage1, as each well is adding ALQ, the current group ALQ needs
1018  // to be communicated to the other processes. But there is no common
1019  // synchronization point that all process will reach in the
1020  // runOptimizeLoop_() in GasLiftSingleWell.cpp.
1021  //
1022  // TODO: Maybe a better solution could be invented by distributing
1023  // wells according to certain parent groups. Then updated group rates
1024  // might not have to be communicated to the other processors.
1025 
1026  // Currently, the best option seems to be to run this part sequentially
1027  // (not in parallel).
1028  //
1029  // TODO: The simplest approach seems to be if a) one process could take
1030  // ownership of all the wells (the union of all the wells in the
1031  // well_container_ of each process) then this process could do the
1032  // optimization, while the other processes could wait for it to
1033  // finish (e.g. comm.barrier()), or alternatively, b) if all
1034  // processes could take ownership of all the wells. Then there
1035  // would be no need for synchronization here..
1036  //
1037  for (int i = 0; i< num_procs; i++) {
1038  int num_rates_to_sync = 0; // communication variable
1039  GLiftSyncGroups groups_to_sync;
1040  if (comm.rank() == i) {
1041  // Run stage1: Optimize single wells while also checking group limits
1042  for (const auto& well : well_container_) {
1043  // NOTE: Only the wells in "group_info" needs to be optimized
1044  if (group_info.hasWell(well->name())) {
1045  gasLiftOptimizationStage1SingleWell(
1046  well.get(), deferred_logger, prod_wells, glift_wells,
1047  group_info, state_map, groups_to_sync
1048  );
1049  }
1050  }
1051  num_rates_to_sync = groups_to_sync.size();
1052  }
1053  // Since "group_info" is not used in stage2, there is no need to
1054  // communicate rates if this is the last iteration...
1055  if (i == (num_procs - 1))
1056  break;
1057  num_rates_to_sync = comm.sum(num_rates_to_sync);
1058  if (num_rates_to_sync > 0) {
1059  std::vector<int> group_indexes;
1060  group_indexes.reserve(num_rates_to_sync);
1061  std::vector<double> group_alq_rates;
1062  group_alq_rates.reserve(num_rates_to_sync);
1063  std::vector<double> group_oil_rates;
1064  group_oil_rates.reserve(num_rates_to_sync);
1065  std::vector<double> group_gas_rates;
1066  group_gas_rates.reserve(num_rates_to_sync);
1067  std::vector<double> group_water_rates;
1068  group_water_rates.reserve(num_rates_to_sync);
1069  if (comm.rank() == i) {
1070  for (auto idx : groups_to_sync) {
1071  auto [oil_rate, gas_rate, water_rate, alq] = group_info.getRates(idx);
1072  group_indexes.push_back(idx);
1073  group_oil_rates.push_back(oil_rate);
1074  group_gas_rates.push_back(gas_rate);
1075  group_water_rates.push_back(water_rate);
1076  group_alq_rates.push_back(alq);
1077  }
1078  } else {
1079  group_indexes.resize(num_rates_to_sync);
1080  group_oil_rates.resize(num_rates_to_sync);
1081  group_gas_rates.resize(num_rates_to_sync);
1082  group_water_rates.resize(num_rates_to_sync);
1083  group_alq_rates.resize(num_rates_to_sync);
1084  }
1085  // TODO: We only need to broadcast to processors with index
1086  // j > i since we do not use the "group_info" in stage 2. In
1087  // this case we should use comm.send() and comm.receive()
1088  // instead of comm.broadcast() to communicate with specific
1089  // processes, and these processes only need to receive the
1090  // data if they are going to check the group rates in stage1
1091  // Another similar idea is to only communicate the rates to
1092  // process j = i + 1
1093 #if HAVE_MPI
1094  EclMpiSerializer ser(comm);
1095  ser.broadcast(i, group_indexes, group_oil_rates,
1096  group_gas_rates, group_water_rates, group_alq_rates);
1097 #endif
1098  if (comm.rank() != i) {
1099  for (int j=0; j<num_rates_to_sync; j++) {
1100  group_info.updateRate(group_indexes[j],
1101  group_oil_rates[j], group_gas_rates[j], group_water_rates[j], group_alq_rates[j]);
1102  }
1103  }
1104  if (this->glift_debug) {
1105  int counter = 0;
1106  if (comm.rank() == i) {
1107  counter = this->wellState().gliftGetDebugCounter();
1108  }
1109  counter = comm.sum(counter);
1110  if (comm.rank() != i) {
1111  this->wellState().gliftSetDebugCounter(counter);
1112  }
1113  }
1114  }
1115  }
1116  }
1117 
1118  // NOTE: this method cannot be const since it passes this->wellState()
1119  // (see below) to the GasLiftSingleWell constructor which accepts WellState
1120  // as a non-const reference..
1121  template<typename TypeTag>
1122  void
1123  BlackoilWellModel<TypeTag>::
1124  gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
1125  DeferredLogger& deferred_logger,
1126  GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
1127  GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
1128  GLiftSyncGroups& sync_groups)
1129  {
1130  const auto& summary_state = ebosSimulator_.vanguard().summaryState();
1131  std::unique_ptr<GasLiftSingleWell> glift
1132  = std::make_unique<GasLiftSingleWell>(
1133  *well, ebosSimulator_, summary_state,
1134  deferred_logger, this->wellState(), this->groupState(),
1135  group_info, sync_groups, this->comm_, this->glift_debug);
1136  auto state = glift->runOptimize(
1137  ebosSimulator_.model().newtonMethod().numIterations());
1138  if (state) {
1139  state_map.insert({well->name(), std::move(state)});
1140  glift_wells.insert({well->name(), std::move(glift)});
1141  return;
1142  }
1143  prod_wells.insert({well->name(), well});
1144  }
1145 
1146 
1147  template<typename TypeTag>
1148  void
1149  BlackoilWellModel<TypeTag>::
1150  initGliftEclWellMap(GLiftEclWells &ecl_well_map)
1151  {
1152  for ( const auto& well: well_container_ ) {
1153  ecl_well_map.try_emplace(
1154  well->name(), &(well->wellEcl()), well->indexOfWell());
1155  }
1156  }
1157 
1158 
1159  template<typename TypeTag>
1160  void
1161  BlackoilWellModel<TypeTag>::
1162  assembleWellEq(const double dt, DeferredLogger& deferred_logger)
1163  {
1164  for (auto& well : well_container_) {
1165  well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1166  }
1167  }
1168 
1169  template<typename TypeTag>
1170  void
1171  BlackoilWellModel<TypeTag>::
1172  apply( BVector& r) const
1173  {
1174  for (auto& well : well_container_) {
1175  well->apply(r);
1176  }
1177  }
1178 
1179 
1180  // Ax = A x - C D^-1 B x
1181  template<typename TypeTag>
1182  void
1183  BlackoilWellModel<TypeTag>::
1184  apply(const BVector& x, BVector& Ax) const
1185  {
1186  for (auto& well : well_container_) {
1187  well->apply(x, Ax);
1188  }
1189  }
1190 
1191  template<typename TypeTag>
1192  void
1193  BlackoilWellModel<TypeTag>::
1194  getWellContributions(WellContributions& wellContribs) const
1195  {
1196  // prepare for StandardWells
1197  wellContribs.setBlockSize(StandardWell<TypeTag>::Indices::numEq, StandardWell<TypeTag>::numStaticWellEq);
1198 
1199  for(unsigned int i = 0; i < well_container_.size(); i++){
1200  auto& well = well_container_[i];
1201  std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1202  if (derived) {
1203  unsigned int numBlocks;
1204  derived->getNumBlocks(numBlocks);
1205  wellContribs.addNumBlocks(numBlocks);
1206  }
1207  }
1208 
1209  // allocate memory for data from StandardWells
1210  wellContribs.alloc();
1211 
1212  for(unsigned int i = 0; i < well_container_.size(); i++){
1213  auto& well = well_container_[i];
1214  // maybe WellInterface could implement addWellContribution()
1215  auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1216  if (derived_std) {
1217  derived_std->addWellContribution(wellContribs);
1218  } else {
1219  auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
1220  if (derived_ms) {
1221  derived_ms->addWellContribution(wellContribs);
1222  } else {
1223  OpmLog::warning("Warning unknown type of well");
1224  }
1225  }
1226  }
1227  }
1228 
1229  // Ax = Ax - alpha * C D^-1 B x
1230  template<typename TypeTag>
1231  void
1232  BlackoilWellModel<TypeTag>::
1233  applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
1234  {
1235  if (this->well_container_.empty()) {
1236  return;
1237  }
1238 
1239  if( scaleAddRes_.size() != Ax.size() ) {
1240  scaleAddRes_.resize( Ax.size() );
1241  }
1242 
1243  scaleAddRes_ = 0.0;
1244  // scaleAddRes_ = - C D^-1 B x
1245  apply( x, scaleAddRes_ );
1246  // Ax = Ax + alpha * scaleAddRes_
1247  Ax.axpy( alpha, scaleAddRes_ );
1248  }
1249 
1250  template<typename TypeTag>
1251  void
1252  BlackoilWellModel<TypeTag>::
1253  addWellContributions(SparseMatrixAdapter& jacobian) const
1254  {
1255  for ( const auto& well: well_container_ ) {
1256  well->addWellContributions(jacobian);
1257  }
1258  }
1259 
1260  template<typename TypeTag>
1261  void
1262  BlackoilWellModel<TypeTag>::
1263  addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const
1264  {
1265  int nw = this->numLocalWellsEnd();
1266  int rdofs = local_num_cells_;
1267  for ( int i = 0; i < nw; i++ ){
1268  int wdof = rdofs + i;
1269  jacobian[wdof][wdof] = 1.0;// better scaling ?
1270  }
1271 
1272  for ( const auto& well : well_container_ ) {
1273  well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
1274  }
1275  }
1276 
1277  template<typename TypeTag>
1278  int
1279  BlackoilWellModel<TypeTag>::
1280  numLocalWellsEnd() const
1281  {
1282  auto w = schedule().getWellsatEnd();
1283  w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end());
1284  return w.size();
1285  }
1286 
1287  template<typename TypeTag>
1288  std::vector<std::vector<int>>
1289  BlackoilWellModel<TypeTag>::
1290  getMaxWellConnections() const
1291  {
1292  std::vector<std::vector<int>> wells;
1293 
1294  auto schedule_wells = schedule().getWellsatEnd();
1295  schedule_wells.erase(std::remove_if(schedule_wells.begin(), schedule_wells.end(), not_on_process_), schedule_wells.end());
1296  wells.reserve(schedule_wells.size());
1297 
1298  // initialize the additional cell connections introduced by wells.
1299  for ( const auto& well : schedule_wells )
1300  {
1301  std::vector<int> compressed_well_perforations;
1302  // All possible completions of the well
1303  const auto& completionSet = well.getConnections();
1304  compressed_well_perforations.reserve(completionSet.size());
1305 
1306  for (const auto& connection: well.getConnections())
1307  {
1308  const int compressed_idx = compressedIndexForInterior(connection.global_index());
1309  if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells.
1310  {
1311  compressed_well_perforations.push_back(compressed_idx);
1312  }
1313  }
1314 
1315  // also include wells with no perforations in case
1316  std::sort(compressed_well_perforations.begin(),
1317  compressed_well_perforations.end());
1318 
1319  wells.push_back(compressed_well_perforations);
1320  }
1321  return wells;
1322  }
1323 
1324  template<typename TypeTag>
1325  void
1326  BlackoilWellModel<TypeTag>::
1327  addWellPressureEquationsStruct(PressureMatrix& jacobian) const
1328  {
1329  int nw = this->numLocalWellsEnd();
1330  int rdofs = local_num_cells_;
1331  for(int i=0; i < nw; i++){
1332  int wdof = rdofs + i;
1333  jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
1334  }
1335  std::vector<std::vector<int>> wellconnections = getMaxWellConnections();
1336  for(int i=0; i < nw; i++){
1337  const auto& perfcells = wellconnections[i];
1338  for(int perfcell : perfcells){
1339  int wdof = rdofs + i;
1340  jacobian.entry(wdof,perfcell) = 0.0;
1341  jacobian.entry(perfcell, wdof) = 0.0;
1342  }
1343  }
1344  }
1345 
1346  template<typename TypeTag>
1347  int
1348  BlackoilWellModel<TypeTag>::
1349  numLocalNonshutWells() const
1350  {
1351  return well_container_.size();
1352  }
1353 
1354  template<typename TypeTag>
1355  void
1356  BlackoilWellModel<TypeTag>::
1357  recoverWellSolutionAndUpdateWellState(const BVector& x)
1358  {
1359 
1360  DeferredLogger local_deferredLogger;
1361  OPM_BEGIN_PARALLEL_TRY_CATCH();
1362  {
1363  for (auto& well : well_container_) {
1364  well->recoverWellSolutionAndUpdateWellState(x, this->wellState(), local_deferredLogger);
1365  }
1366 
1367  }
1368  OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1369  "recoverWellSolutionAndUpdateWellState() failed: ",
1370  terminal_output_, ebosSimulator_.vanguard().grid().comm());
1371 
1372  }
1373 
1374 
1375 
1376 
1377  template<typename TypeTag>
1378  void
1379  BlackoilWellModel<TypeTag>::
1380  initPrimaryVariablesEvaluation() const
1381  {
1382  for (auto& well : well_container_) {
1383  well->initPrimaryVariablesEvaluation();
1384  }
1385  }
1386 
1387 
1388 
1389 
1390 
1391 
1392  template<typename TypeTag>
1393  ConvergenceReport
1394  BlackoilWellModel<TypeTag>::
1395  getWellConvergence(const std::vector<Scalar>& B_avg, bool checkWellGroupControls) const
1396  {
1397 
1398  DeferredLogger local_deferredLogger;
1399  // Get global (from all processes) convergence report.
1400  ConvergenceReport local_report;
1401  const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1402  for (const auto& well : well_container_) {
1403  if (well->isOperableAndSolvable() || well->wellIsStopped()) {
1404  local_report += well->getWellConvergence(this->wellState(), B_avg, local_deferredLogger, iterationIdx > param_.strict_outer_iter_wells_ );
1405  } else {
1406  ConvergenceReport report;
1407  using CR = ConvergenceReport;
1408  report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
1409  local_report += report;
1410  }
1411  }
1412 
1413  const Opm::Parallel::Communication comm = grid().comm();
1414  DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
1415  ConvergenceReport report = gatherConvergenceReport(local_report, comm);
1416 
1417  // the well_group_control_changed info is already communicated
1418  if (checkWellGroupControls) {
1419  report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed);
1420  }
1421 
1422  if (terminal_output_) {
1423  global_deferredLogger.logMessages();
1424  }
1425  // Log debug messages for NaN or too large residuals.
1426  if (terminal_output_) {
1427  for (const auto& f : report.wellFailures()) {
1428  if (f.severity() == ConvergenceReport::Severity::NotANumber) {
1429  OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1430  } else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
1431  OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1432  }
1433  }
1434  }
1435  return report;
1436  }
1437 
1438 
1439 
1440 
1441 
1442  template<typename TypeTag>
1443  void
1445  calculateExplicitQuantities(DeferredLogger& deferred_logger) const
1446  {
1447  // TODO: checking isOperableAndSolvable() ?
1448  for (auto& well : well_container_) {
1449  well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), deferred_logger);
1450  }
1451  }
1452 
1453 
1454 
1455 
1456 
1457  template<typename TypeTag>
1458  bool
1460  shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const
1461  {
1462  const auto& balance = schedule()[reportStepIdx].network_balance();
1463  if (balance.mode() == Network::Balance::CalcMode::TimeStepStart) {
1464  return iterationIdx == 0;
1465  } else if (balance.mode() == Network::Balance::CalcMode::NUPCOL) {
1466  const int nupcol = schedule()[reportStepIdx].nupcol();
1467  return iterationIdx < nupcol;
1468  } else {
1469  // We do not support any other rebalancing modes,
1470  // i.e. TimeInterval based rebalancing is not available.
1471  // This should be warned about elsewhere, so we choose to
1472  // avoid spamming with a warning here.
1473  return false;
1474  }
1475  }
1476 
1477 
1478 
1479 
1480 
1481  template<typename TypeTag>
1482  std::tuple<bool, bool, double>
1483  BlackoilWellModel<TypeTag>::
1484  updateWellControls(DeferredLogger& deferred_logger)
1485  {
1486  // Even if there are no wells active locally, we cannot
1487  // return as the DeferredLogger uses global communication.
1488  // For no well active globally we simply return.
1489  if( !wellsActive() ) return { false, false, 0.0 };
1490 
1491  const int episodeIdx = ebosSimulator_.episodeIndex();
1492  const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1493  const auto& comm = ebosSimulator_.vanguard().grid().comm();
1494  updateAndCommunicateGroupData(episodeIdx, iterationIdx);
1495 
1496  const auto [local_network_changed, local_network_imbalance]
1497  = shouldBalanceNetwork(episodeIdx, iterationIdx) ?
1498  updateNetworkPressures(episodeIdx) : std::make_pair(false, 0.0);
1499  const bool network_changed = comm.sum(local_network_changed);
1500  const double network_imbalance = comm.max(local_network_imbalance);
1501 
1502  bool changed_well_group = false;
1503  // Check group individual constraints.
1504  const int nupcol = schedule()[episodeIdx].nupcol();
1505  // don't switch group control when iterationIdx > nupcol
1506  // to avoid oscilations between group controls
1507  if (iterationIdx <= nupcol) {
1508  const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
1509  changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
1510  }
1511  // Check wells' group constraints and communicate.
1512  bool changed_well_to_group = false;
1513  for (const auto& well : well_container_) {
1514  const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
1515  const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1516  if (changed_well) {
1517  changed_well_to_group = changed_well || changed_well_to_group;
1518  }
1519  }
1520 
1521  changed_well_to_group = comm.sum(changed_well_to_group);
1522  if (changed_well_to_group) {
1523  updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1524  changed_well_group = true;
1525  }
1526 
1527  // Check individual well constraints and communicate.
1528  bool changed_well_individual = false;
1529  for (const auto& well : well_container_) {
1530  const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
1531  const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1532  if (changed_well) {
1533  changed_well_individual = changed_well || changed_well_individual;
1534  }
1535  }
1536  changed_well_individual = comm.sum(changed_well_individual);
1537  if (changed_well_individual) {
1538  updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1539  changed_well_group = true;
1540  }
1541 
1542  // update wsolvent fraction for REIN wells
1543  const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
1544  updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState());
1545 
1546  return { changed_well_group, network_changed, network_imbalance };
1547  }
1548 
1549 
1550  template<typename TypeTag>
1551  void
1552  BlackoilWellModel<TypeTag>::
1553  updateAndCommunicate(const int reportStepIdx,
1554  const int iterationIdx,
1555  DeferredLogger& deferred_logger)
1556  {
1557  updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
1558  // if a well or group change control it affects all wells that are under the same group
1559  for (const auto& well : well_container_) {
1560  well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
1561  }
1562  updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
1563  }
1564 
1565  template<typename TypeTag>
1566  bool
1567  BlackoilWellModel<TypeTag>::
1568  updateGroupControls(const Group& group,
1569  DeferredLogger& deferred_logger,
1570  const int reportStepIdx,
1571  const int iterationIdx)
1572  {
1573  bool changed = false;
1574  bool changed_hc = checkGroupHigherConstraints( group, deferred_logger, reportStepIdx);
1575  if (changed_hc) {
1576  changed = true;
1577  updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1578  }
1579  bool changed_individual = updateGroupIndividualControl( group, deferred_logger, reportStepIdx);
1580  if (changed_individual) {
1581  changed = true;
1582  updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1583  }
1584  // call recursively down the group hierarchy
1585  for (const std::string& groupName : group.groups()) {
1586  bool changed_this = updateGroupControls( schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx,iterationIdx);
1587  changed = changed || changed_this;
1588  }
1589  return changed;
1590  }
1591 
1592  template<typename TypeTag>
1593  void
1595  updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
1596  {
1597  DeferredLogger local_deferredLogger;
1598  for (const auto& well : well_container_) {
1599  const auto& wname = well->name();
1600  const auto wasClosed = wellTestState.well_is_closed(wname);
1601  well->checkWellOperability(ebosSimulator_, this->wellState(), local_deferredLogger);
1602  well->updateWellTestState(this->wellState().well(wname), simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
1603 
1604  if (!wasClosed && wellTestState.well_is_closed(wname)) {
1605  this->closed_this_step_.insert(wname);
1606  }
1607  }
1608 
1609  const Opm::Parallel::Communication comm = grid().comm();
1610  DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
1611  if (terminal_output_) {
1612  global_deferredLogger.logMessages();
1613  }
1614  }
1615 
1616 
1617  template<typename TypeTag>
1618  void
1619  BlackoilWellModel<TypeTag>::computePotentials(const std::size_t widx,
1620  const WellState& well_state_copy,
1621  std::string& exc_msg,
1622  ExceptionType::ExcEnum& exc_type,
1623  DeferredLogger& deferred_logger)
1624  {
1625  const int np = numPhases();
1626  std::vector<double> potentials;
1627  const auto& well= well_container_[widx];
1628  try {
1629  well->computeWellPotentials(ebosSimulator_, well_state_copy, potentials, deferred_logger);
1630  }
1631  // catch all possible exception and store type and message.
1632  OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
1633  // Store it in the well state
1634  // potentials is resized and set to zero in the beginning of well->ComputeWellPotentials
1635  // and updated only if sucessfull. i.e. the potentials are zero for exceptions
1636  auto& ws = this->wellState().well(well->indexOfWell());
1637  for (int p = 0; p < np; ++p) {
1638  // make sure the potentials are positive
1639  ws.well_potentials[p] = std::max(0.0, potentials[p]);
1640  }
1641  }
1642 
1643 
1644 
1645  template <typename TypeTag>
1646  void
1647  BlackoilWellModel<TypeTag>::
1648  calculateProductivityIndexValues(DeferredLogger& deferred_logger)
1649  {
1650  for (const auto& wellPtr : this->well_container_) {
1651  this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
1652  }
1653  }
1654 
1655 
1656 
1657 
1658 
1659  template <typename TypeTag>
1660  void
1661  BlackoilWellModel<TypeTag>::
1662  calculateProductivityIndexValuesShutWells(const int reportStepIdx,
1663  DeferredLogger& deferred_logger)
1664  {
1665  // For the purpose of computing PI/II values, it is sufficient to
1666  // construct StandardWell instances only. We don't need to form
1667  // well objects that honour the 'isMultisegment()' flag of the
1668  // corresponding "this->wells_ecl_[shutWell]".
1669 
1670  for (const auto& shutWell : this->local_shut_wells_) {
1671  if (this->wells_ecl_[shutWell].getConnections().empty()) {
1672  // No connections in this well. Nothing to do.
1673  continue;
1674  }
1675 
1676  auto wellPtr = this->template createTypedWellPointer
1677  <StandardWell<TypeTag>>(shutWell, reportStepIdx);
1678 
1679  wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
1680  this->local_num_cells_, this->B_avg_, true);
1681 
1682  this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
1683  }
1684  }
1685 
1686 
1687 
1688 
1689 
1690  template <typename TypeTag>
1691  void
1692  BlackoilWellModel<TypeTag>::
1693  calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
1694  DeferredLogger& deferred_logger)
1695  {
1696  wellPtr->updateProductivityIndex(this->ebosSimulator_,
1697  this->prod_index_calc_[wellPtr->indexOfWell()],
1698  this->wellState(),
1699  deferred_logger);
1700  }
1701 
1702 
1703 
1704 
1705 
1706  template<typename TypeTag>
1707  void
1708  BlackoilWellModel<TypeTag>::
1709  prepareTimeStep(DeferredLogger& deferred_logger)
1710  {
1711  for (const auto& well : well_container_) {
1712  auto& events = this->wellState().well(well->indexOfWell()).events;
1713  if (events.hasEvent(WellState::event_mask)) {
1714  well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
1715  well->updatePrimaryVariables(this->wellState(), deferred_logger);
1716  well->initPrimaryVariablesEvaluation();
1717  // There is no new well control change input within a report step,
1718  // so next time step, the well does not consider to have effective events anymore.
1719  events.clearEvent(WellState::event_mask);
1720  }
1721  // solve the well equation initially to improve the initial solution of the well model
1722  if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
1723  try {
1724  well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), deferred_logger);
1725  } catch (const std::exception& e) {
1726  const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the privious rates";
1727  deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
1728  }
1729  }
1730  }
1731  updatePrimaryVariables(deferred_logger);
1732  }
1733 
1734  template<typename TypeTag>
1735  void
1736  BlackoilWellModel<TypeTag>::
1737  updateAverageFormationFactor()
1738  {
1739  std::vector< Scalar > B_avg(numComponents(), Scalar() );
1740  const auto& grid = ebosSimulator_.vanguard().grid();
1741  const auto& gridView = grid.leafGridView();
1742  ElementContext elemCtx(ebosSimulator_);
1743 
1744  OPM_BEGIN_PARALLEL_TRY_CATCH();
1745  for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1746  elemCtx.updatePrimaryStencil(elem);
1747  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1748 
1749  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
1750  const auto& fs = intQuants.fluidState();
1751 
1752  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1753  {
1754  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1755  continue;
1756  }
1757 
1758  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1759  auto& B = B_avg[ compIdx ];
1760 
1761  B += 1 / fs.invB(phaseIdx).value();
1762  }
1763  if constexpr (has_solvent_) {
1764  auto& B = B_avg[solventSaturationIdx];
1765  B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
1766  }
1767  }
1768  OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updateAverageFormationFactor() failed: ", grid.comm())
1769 
1770  // compute global average
1771  grid.comm().sum(B_avg.data(), B_avg.size());
1772  for (auto& bval : B_avg)
1773  {
1774  bval /= global_num_cells_;
1775  }
1776  B_avg_ = B_avg;
1777  }
1778 
1779 
1780 
1781 
1782 
1783  template<typename TypeTag>
1784  void
1785  BlackoilWellModel<TypeTag>::
1786  updatePrimaryVariables(DeferredLogger& deferred_logger)
1787  {
1788  for (const auto& well : well_container_) {
1789  well->updatePrimaryVariables(this->wellState(), deferred_logger);
1790  }
1791  }
1792 
1793  template<typename TypeTag>
1794  void
1795  BlackoilWellModel<TypeTag>::extractLegacyCellPvtRegionIndex_()
1796  {
1797  const auto& grid = ebosSimulator_.vanguard().grid();
1798  const auto& eclProblem = ebosSimulator_.problem();
1799  const unsigned numCells = grid.size(/*codim=*/0);
1800 
1801  pvt_region_idx_.resize(numCells);
1802  for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
1803  pvt_region_idx_[cellIdx] =
1804  eclProblem.pvtRegionIndex(cellIdx);
1805  }
1806  }
1807 
1808  // The number of components in the model.
1809  template<typename TypeTag>
1810  int
1811  BlackoilWellModel<TypeTag>::numComponents() const
1812  {
1813  // The numComponents here does not reflect the actual number of the components in the system.
1814  // It more or less reflects the number of mass conservation equations for the well equations.
1815  // For example, in the current formulation, we do not have the polymer conservation equation
1816  // in the well equations. As a result, for an oil-water-polymer system, this function will return 2.
1817  // In some way, it makes this function appear to be confusing from its name, and we need
1818  // to revisit/revise this function again when extending the variants of system that flow can simulate.
1819  if (numPhases() < 3) {
1820  return numPhases();
1821  }
1822  int numComp = FluidSystem::numComponents;
1823  if constexpr (has_solvent_) {
1824  numComp ++;
1825  }
1826 
1827  return numComp;
1828  }
1829 
1830  template<typename TypeTag>
1831  void
1832  BlackoilWellModel<TypeTag>::extractLegacyDepth_()
1833  {
1834  const auto& eclProblem = ebosSimulator_.problem();
1835  depth_.resize(local_num_cells_);
1836  for (unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
1837  depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
1838  }
1839  }
1840 
1841  template<typename TypeTag>
1842  void
1843  BlackoilWellModel<TypeTag>::
1844  updatePerforationIntensiveQuantities()
1845  {
1846  ElementContext elemCtx(ebosSimulator_);
1847  const auto& gridView = ebosSimulator_.gridView();
1848 
1849  OPM_BEGIN_PARALLEL_TRY_CATCH();
1850  for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1851  elemCtx.updatePrimaryStencil(elem);
1852  int elemIdx = elemCtx.globalSpaceIndex(0, 0);
1853 
1854  if (!is_cell_perforated_[elemIdx]) {
1855  continue;
1856  }
1857  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1858  }
1859  OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updatePerforationIntensiveQuantities() failed: ", ebosSimulator_.vanguard().grid().comm());
1860  }
1861 
1862 
1863  template<typename TypeTag>
1864  typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1865  BlackoilWellModel<TypeTag>::
1866  getWell(const std::string& well_name) const
1867  {
1868  // finding the iterator of the well in wells_ecl
1869  auto well = std::find_if(well_container_.begin(),
1870  well_container_.end(),
1871  [&well_name](const WellInterfacePtr& elem)->bool {
1872  return elem->name() == well_name;
1873  });
1874 
1875  assert(well != well_container_.end());
1876 
1877  return *well;
1878  }
1879 
1880  template<typename TypeTag>
1881  bool
1882  BlackoilWellModel<TypeTag>::
1883  hasWell(const std::string& well_name) const
1884  {
1885  return std::any_of(well_container_.begin(), well_container_.end(),
1886  [&well_name](const WellInterfacePtr& elem) -> bool
1887  {
1888  return elem->name() == well_name;
1889  });
1890  }
1891 
1892 
1893 
1894 
1895  template <typename TypeTag>
1896  int
1897  BlackoilWellModel<TypeTag>::
1898  reportStepIndex() const
1899  {
1900  return std::max(this->ebosSimulator_.episodeIndex(), 0);
1901  }
1902 
1903 
1904 
1905 
1906 
1907  template<typename TypeTag>
1908  void
1909  BlackoilWellModel<TypeTag>::
1910  calcRates(const int fipnum,
1911  const int pvtreg,
1912  std::vector<double>& resv_coeff)
1913  {
1914  rateConverter_->calcCoeff(fipnum, pvtreg, resv_coeff);
1915  }
1916 
1917  template<typename TypeTag>
1918  void
1919  BlackoilWellModel<TypeTag>::
1920  calcInjRates(const int fipnum,
1921  const int pvtreg,
1922  std::vector<double>& resv_coeff)
1923  {
1924  rateConverter_->calcInjCoeff(fipnum, pvtreg, resv_coeff);
1925  }
1926 
1927 
1928  template <typename TypeTag>
1929  void
1930  BlackoilWellModel<TypeTag>::
1931  computeWellTemperature()
1932  {
1933  if (!has_energy_)
1934  return;
1935 
1936  int np = numPhases();
1937  double cellInternalEnergy;
1938  double cellBinv;
1939  double cellDensity;
1940  double perfPhaseRate;
1941  const int nw = numLocalWells();
1942  for (auto wellID = 0*nw; wellID < nw; ++wellID) {
1943  const Well& well = wells_ecl_[wellID];
1944  if (well.isInjector())
1945  continue;
1946 
1947  int connpos = 0;
1948  for (int i = 0; i < wellID; ++i) {
1949  connpos += well_perf_data_[i].size();
1950  }
1951  connpos *= np;
1952  std::array<double,2> weighted{0.0,0.0};
1953  auto& [weighted_temperature, total_weight] = weighted;
1954 
1955  auto& well_info = local_parallel_well_info_[wellID].get();
1956  const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size());
1957  auto& ws = this->wellState().well(wellID);
1958  auto& perf_data = ws.perf_data;
1959  auto& perf_phase_rate = perf_data.phase_rates;
1960 
1961  for (int perf = 0; perf < num_perf_this_well; ++perf) {
1962  const int cell_idx = well_perf_data_[wellID][perf].cell_index;
1963  const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1964  const auto& fs = intQuants.fluidState();
1965 
1966  // we on only have one temperature pr cell any phaseIdx will do
1967  double cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
1968 
1969  double weight_factor = 0.0;
1970  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1971  {
1972  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1973  continue;
1974  }
1975  cellInternalEnergy = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
1976  cellBinv = fs.invB(phaseIdx).value();
1977  cellDensity = fs.density(phaseIdx).value();
1978  perfPhaseRate = perf_phase_rate[ perf*np + phaseIdx ];
1979  weight_factor += cellDensity * perfPhaseRate/cellBinv * cellInternalEnergy/cellTemperatures;
1980  }
1981  total_weight += weight_factor;
1982  weighted_temperature += weight_factor * cellTemperatures;
1983  }
1984  well_info.communication().sum(weighted.data(), 2);
1985  this->wellState().well(wellID).temperature = weighted_temperature/total_weight;
1986  }
1987  }
1988 
1989 
1990 
1991  template <typename TypeTag>
1992  void
1993  BlackoilWellModel<TypeTag>::
1994  assignWellTracerRates(data::Wells& wsrpt) const
1995  {
1996  const auto & wellTracerRates = ebosSimulator_.problem().tracerModel().getWellTracerRates();
1997 
1998  if (wellTracerRates.empty())
1999  return; // no tracers
2000 
2001  for (const auto& wTR : wellTracerRates) {
2002  std::string wellName = wTR.first.first;
2003  auto xwPos = wsrpt.find(wellName);
2004  if (xwPos == wsrpt.end()) { // No well results.
2005  continue;
2006  }
2007  std::string tracerName = wTR.first.second;
2008  double rate = wTR.second;
2009  xwPos->second.rates.set(data::Rates::opt::tracer, rate, tracerName);
2010  }
2011  }
2012 
2013 
2014 
2015 
2016 
2017 } // namespace Opm
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:94
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Calculating the explict quantities used in the well calculation.
Definition: BlackoilWellModel_impl.hpp:1445
Definition: DeferredLogger.hpp:57
void logMessages()
Log all messages to the OpmLog backends, and clear the message container.
Definition: DeferredLogger.cpp:85
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Opm::DeferredLogger gatherDeferredLogger(const Opm::DeferredLogger &local_deferredlogger, Opm::Parallel::Communication)
Create a global log combining local logs.
Definition: gatherDeferredLogger.cpp:168
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition: phaseUsageFromDeck.cpp:145