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