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