My Project
Loading...
Searching...
No Matches
BlackoilWellModel_impl.hpp
1/*
2 Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics.
3 Copyright 2016 - 2018 Equinor ASA.
4 Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
5 Copyright 2016 - 2018 Norce AS
6
7 This file is part of the Open Porous Media project (OPM).
8
9 OPM is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 OPM is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with OPM. If not, see <http://www.gnu.org/licenses/>.
21*/
22
23#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24#include <opm/core/props/phaseUsageFromDeck.hpp>
25#include <opm/grid/utility/cartesianToCompressed.hpp>
26
27#include <opm/input/eclipse/Units/UnitSystem.hpp>
28#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
29#include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
30#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
31
32#include <opm/simulators/wells/BlackoilWellModelConstraints.hpp>
33#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
34#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
35#include <opm/simulators/wells/VFPProperties.hpp>
36#include <opm/simulators/utils/MPIPacker.hpp>
37#include <opm/simulators/linalg/bda/WellContributions.hpp>
38
39#if HAVE_MPI
40#include <ebos/eclmpiserializer.hh>
41#endif
42
43#include <algorithm>
44#include <iomanip>
45#include <utility>
46
47#include <fmt/format.h>
48
49namespace Opm {
50 template<typename TypeTag>
51 BlackoilWellModel<TypeTag>::
52 BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage)
53 : BlackoilWellModelGeneric(ebosSimulator.vanguard().schedule(),
54 ebosSimulator.vanguard().summaryState(),
55 ebosSimulator.vanguard().eclState(),
56 phase_usage,
57 ebosSimulator.gridView().comm())
58 , ebosSimulator_(ebosSimulator)
59 {
60 terminal_output_ = ((ebosSimulator.gridView().comm().rank() == 0) &&
61 EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput));
62
63 local_num_cells_ = ebosSimulator_.gridView().size(0);
64
65 // Number of cells the global grid view
66 global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();
67
68 {
69 auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
70
71 this->parallel_well_info_.reserve(parallel_wells.size());
72 for( const auto& name_bool : parallel_wells) {
73 this->parallel_well_info_.emplace_back(name_bool, grid().comm());
74 }
75 }
76
77 this->alternative_well_rate_init_ =
78 EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
79
80 this->wbpCalculationService_
81 .localCellIndex([this](const std::size_t globalIndex)
82 { return this->compressedIndexForInterior(globalIndex); })
83 .evalCellSource([this](const int localCell,
84 PAvgDynamicSourceData::SourceDataSpan<double> sourceTerms)
85 {
86 using Item = PAvgDynamicSourceData::SourceDataSpan<double>::Item;
87
88 const auto* intQuants = this->ebosSimulator_.model()
89 .cachedIntensiveQuantities(localCell, /*timeIndex = */0);
90 const auto& fs = intQuants->fluidState();
91
92 sourceTerms.set(Item::PoreVol, intQuants->porosity().value() *
93 this->ebosSimulator_.model().dofTotalVolume(localCell));
94
95 constexpr auto io = FluidSystem::oilPhaseIdx;
96 constexpr auto ig = FluidSystem::gasPhaseIdx;
97 constexpr auto iw = FluidSystem::waterPhaseIdx;
98
99 // Ideally, these would be 'constexpr'.
100 const auto haveOil = FluidSystem::phaseIsActive(io);
101 const auto haveGas = FluidSystem::phaseIsActive(ig);
102 const auto haveWat = FluidSystem::phaseIsActive(iw);
103
104 auto weightedPhaseDensity = [&fs](const auto ip)
105 {
106 return fs.saturation(ip).value() * fs.density(ip).value();
107 };
108
109 if (haveOil) { sourceTerms.set(Item::Pressure, fs.pressure(io).value()); }
110 else if (haveGas) { sourceTerms.set(Item::Pressure, fs.pressure(ig).value()); }
111 else { sourceTerms.set(Item::Pressure, fs.pressure(iw).value()); }
112
113 // Strictly speaking, assumes SUM(s[p]) == 1.
114 auto rho = 0.0;
115 if (haveOil) { rho += weightedPhaseDensity(io); }
116 if (haveGas) { rho += weightedPhaseDensity(ig); }
117 if (haveWat) { rho += weightedPhaseDensity(iw); }
118
119 sourceTerms.set(Item::MixtureDensity, rho);
120 });
121 }
122
123 template<typename TypeTag>
124 BlackoilWellModel<TypeTag>::
125 BlackoilWellModel(Simulator& ebosSimulator) :
126 BlackoilWellModel(ebosSimulator, phaseUsageFromDeck(ebosSimulator.vanguard().eclState()))
127 {}
128
129
130 template<typename TypeTag>
131 void
132 BlackoilWellModel<TypeTag>::
133 init()
134 {
135 extractLegacyCellPvtRegionIndex_();
136 extractLegacyDepth_();
137
138 gravity_ = ebosSimulator_.problem().gravity()[2];
139
140 initial_step_ = true;
141
142 // add the eWoms auxiliary module for the wells to the list
143 ebosSimulator_.model().addAuxiliaryModule(this);
144
145 is_cell_perforated_.resize(local_num_cells_, false);
146 }
147
148
149 template<typename TypeTag>
150 void
151 BlackoilWellModel<TypeTag>::
152 initWellContainer(const int reportStepIdx)
153 {
154 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
155 + ScheduleEvents::NEW_WELL;
156 const auto& events = schedule()[reportStepIdx].wellgroup_events();
157 for (auto& wellPtr : this->well_container_) {
158 const bool well_opened_this_step = report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask);
159 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
160 this->local_num_cells_, this->B_avg_, well_opened_this_step);
161 }
162 }
163
164
165 template<typename TypeTag>
166 void
167 BlackoilWellModel<TypeTag>::
168 addNeighbors(std::vector<NeighborSet>& neighbors) const
169 {
170 if (!param_.matrix_add_well_contributions_) {
171 return;
172 }
173
174 // Create cartesian to compressed mapping
175 const auto& schedule_wells = schedule().getWellsatEnd();
176
177 // initialize the additional cell connections introduced by wells.
178 for (const auto& well : schedule_wells)
179 {
180 std::vector<int> wellCells = this->getCellsForConnections(well);
181 for (int cellIdx : wellCells) {
182 neighbors[cellIdx].insert(wellCells.begin(),
183 wellCells.end());
184 }
185 }
186 }
187
188
189 template<typename TypeTag>
190 void
191 BlackoilWellModel<TypeTag>::
192 linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
193 {
194 OPM_BEGIN_PARALLEL_TRY_CATCH();
195 for (const auto& well: well_container_) {
196 // Modifiy the Jacobian with explicit Schur complement
197 // contributions if requested.
198 if (param_.matrix_add_well_contributions_) {
199 well->addWellContributions(jacobian);
200 }
201 // Apply as Schur complement the well residual to reservoir residuals:
202 // r = r - duneC_^T * invDuneD_ * resWell_
203 well->apply(res);
204 }
205 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ",
206 ebosSimulator_.gridView().comm());
207 }
208
209
210 template<typename TypeTag>
211 void
212 BlackoilWellModel<TypeTag>::
213 linearizeDomain(const Domain& domain, SparseMatrixAdapter& jacobian, GlobalEqVector& res)
214 {
215 // Note: no point in trying to do a parallel gathering
216 // try/catch here, as this function is not called in
217 // parallel but for each individual domain of each rank.
218 for (const auto& well: well_container_) {
219 if (well_domain_.at(well->name()) == domain.index) {
220 // Modifiy the Jacobian with explicit Schur complement
221 // contributions if requested.
222 if (param_.matrix_add_well_contributions_) {
223 well->addWellContributions(jacobian);
224 }
225 // Apply as Schur complement the well residual to reservoir residuals:
226 // r = r - duneC_^T * invDuneD_ * resWell_
227 well->apply(res);
228 }
229 }
230 }
231
232
233 template<typename TypeTag>
234 void
235 BlackoilWellModel<TypeTag>::
236 beginReportStep(const int timeStepIdx)
237 {
238 DeferredLogger local_deferredLogger{};
239
240 this->report_step_starts_ = true;
241
242 {
243 // WELPI scaling runs at start of report step.
244 const auto enableWellPIScaling = true;
245 this->initializeLocalWellStructure(timeStepIdx, enableWellPIScaling);
246 }
247
248 this->initializeGroupStructure(timeStepIdx);
249
250 const auto& comm = this->ebosSimulator_.vanguard().grid().comm();
251
252 OPM_BEGIN_PARALLEL_TRY_CATCH()
253 {
254 // Create facility for calculating reservoir voidage volumes for
255 // purpose of RESV controls.
256 this->rateConverter_ = std::make_unique<RateConverterType>
257 (this->phase_usage_, std::vector<int>(this->local_num_cells_, 0));
258 this->rateConverter_->template defineState<ElementContext>(this->ebosSimulator_);
259
260 // Update VFP properties.
261 {
262 const auto& sched_state = this->schedule()[timeStepIdx];
263
264 this->vfp_properties_ = std::make_unique<VFPProperties>
265 (sched_state.vfpinj(), sched_state.vfpprod(), this->wellState());
266 }
267 }
268 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
269 "beginReportStep() failed: ",
270 this->terminal_output_, comm)
271
272 // Store the current well and group states in order to recover in
273 // the case of failed iterations
274 this->commitWGState();
275 }
276
277
278
279
280
281 template <typename TypeTag>
282 void
285 const bool enableWellPIScaling)
286 {
288
289 const auto& comm = this->ebosSimulator_.vanguard().grid().comm();
290
291 // Wells_ecl_ holds this rank's wells, both open and stopped/shut.
292 this->wells_ecl_ = this->getLocalWells(reportStepIdx);
293 this->local_parallel_well_info_ =
294 this->createLocalParallelWellInfo(this->wells_ecl_);
295
296 // At least initializeWellState() might be throw an exception in
297 // UniformTabulated2DFunction. Playing it safe by extending the
298 // scope a bit.
299 OPM_BEGIN_PARALLEL_TRY_CATCH()
300 {
301 this->initializeWellPerfData();
302 this->initializeWellState(reportStepIdx);
303 this->initializeWBPCalculationService();
304
305 if (this->param_.use_multisegment_well_ && this->anyMSWellOpenLocal()) {
306 this->wellState().initWellStateMSWell(this->wells_ecl_, &this->prevWellState());
307 }
308
309 this->initializeWellProdIndCalculators();
310
311 if (enableWellPIScaling && this->schedule()[reportStepIdx].events()
312 .hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX))
313 {
314 this->runWellPIScaling(reportStepIdx, local_deferredLogger);
315 }
316 }
317 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
318 "Failed to initialize local well structure: ",
319 this->terminal_output_, comm)
320 }
321
322
323
324
325
326 template <typename TypeTag>
327 void
330 {
332
333 const auto& comm = this->ebosSimulator_.vanguard().grid().comm();
334
335 OPM_BEGIN_PARALLEL_TRY_CATCH()
336 {
337 const auto& fieldGroup =
338 this->schedule().getGroup("FIELD", reportStepIdx);
339
340 WellGroupHelpers::setCmodeGroup(fieldGroup,
341 this->schedule(),
342 this->summaryState(),
344 this->groupState());
345
346 // Define per region average pressure calculators for use by
347 // pressure maintenance groups (GPMAINT keyword).
348 if (this->schedule()[reportStepIdx].has_gpmaint()) {
349 WellGroupHelpers::setRegionAveragePressureCalculator
350 (fieldGroup,
351 this->schedule(),
353 this->eclState_.fieldProps(),
354 this->phase_usage_,
355 this->regionalAveragePressureCalculator_);
356 }
357 }
358 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
359 "Failed to initialize group structure: ",
360 this->terminal_output_, comm)
361 }
362
363
364
365
366
367 // called at the beginning of a time step
368 template<typename TypeTag>
369 void
372 {
373 OPM_TIMEBLOCK(beginTimeStep);
374 updateAverageFormationFactor();
376 switched_prod_groups_.clear();
377 switched_inj_groups_.clear();
378
379 this->resetWGState();
380 const int reportStepIdx = ebosSimulator_.episodeIndex();
381 updateAndCommunicateGroupData(reportStepIdx,
382 ebosSimulator_.model().newtonMethod().numIterations());
383 this->wellState().gliftTimeStepInit();
384 const double simulationTime = ebosSimulator_.time();
385 OPM_BEGIN_PARALLEL_TRY_CATCH();
386 {
387 // test wells
389
390 // create the well container
391 createWellContainer(reportStepIdx);
392
393 // Wells are active if they are active wells on at least one process.
394 const Grid& grid = ebosSimulator_.vanguard().grid();
395 wells_active_ = !this->well_container_.empty();
396 wells_active_ = grid.comm().max(wells_active_);
397
398 // do the initialization for all the wells
399 // TODO: to see whether we can postpone of the intialization of the well containers to
400 // optimize the usage of the following several member variables
401 this->initWellContainer(reportStepIdx);
402
403 // update the updated cell flag
404 std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
405 for (auto& well : well_container_) {
406 well->updatePerforatedCell(is_cell_perforated_);
407 }
408
409 // calculate the efficiency factors for each well
410 calculateEfficiencyFactors(reportStepIdx);
411
412 if constexpr (has_polymer_)
413 {
414 if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
415 setRepRadiusPerfLength();
416 }
417 }
418
419 }
420 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "beginTimeStep() failed: ",
421 terminal_output_, ebosSimulator_.vanguard().grid().comm());
422
423 for (auto& well : well_container_) {
424 well->setVFPProperties(vfp_properties_.get());
425 well->setGuideRate(&guideRate_);
426 }
427
428 this->updateInjFCMult(local_deferredLogger);
429
430 // Close completions due to economic reasons
431 for (auto& well : well_container_) {
432 well->closeCompletions(wellTestState());
433 }
434
435 // we need the inj_multiplier from the previous time step
436 this->initInjMult();
437
438 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
439 if (alternative_well_rate_init_) {
440 // Update the well rates of well_state_, if only single-phase rates, to
441 // have proper multi-phase rates proportional to rates at bhp zero.
442 // This is done only for producers, as injectors will only have a single
443 // nonzero phase anyway.
444 for (auto& well : well_container_) {
445 const bool zero_target = well->stopppedOrZeroRateTarget(summaryState, this->wellState());
446 if (well->isProducer() && !zero_target) {
447 well->updateWellStateRates(ebosSimulator_, this->wellState(), local_deferredLogger);
448 }
449 }
450 }
451
452 for (auto& well : well_container_) {
453 if (well->isVFPActive(local_deferredLogger)){
454 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
455 }
456 }
457
458 // calculate the well potentials
459 try {
460 updateWellPotentials(reportStepIdx,
461 /*onlyAfterEvent*/true,
462 ebosSimulator_.vanguard().summaryConfig(),
463 local_deferredLogger);
464 } catch ( std::runtime_error& e ) {
465 const std::string msg = "A zero well potential is returned for output purposes. ";
466 local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
467 }
468
469 //update guide rates
470 const auto& comm = ebosSimulator_.vanguard().grid().comm();
471 std::vector<double> pot(numPhases(), 0.0);
472 const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
473 WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
474 this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
475 std::string exc_msg;
476 auto exc_type = ExceptionType::NONE;
477 // update gpmaint targets
478 if (schedule_[reportStepIdx].has_gpmaint()) {
479 for (auto& calculator : regionalAveragePressureCalculator_) {
480 calculator.second->template defineState<ElementContext>(ebosSimulator_);
481 }
482 const double dt = ebosSimulator_.timeStepSize();
483 WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup,
484 schedule_, regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState());
485 }
486 try {
487 // Compute initial well solution for new wells and injectors that change injection type i.e. WAG.
488 for (auto& well : well_container_) {
489 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
490 + ScheduleEvents::INJECTION_TYPE_CHANGED
491 + ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
492 + ScheduleEvents::NEW_WELL;
493
494 const auto& events = schedule()[reportStepIdx].wellgroup_events();
495 const bool event = report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
496 const bool dyn_status_change = this->wellState().well(well->name()).status
497 != this->prevWellState().well(well->name()).status;
498
499 if (event || dyn_status_change) {
500 try {
501 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), local_deferredLogger);
502 well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), local_deferredLogger);
503 well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), local_deferredLogger);
504 } catch (const std::exception& e) {
505 const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates";
506 local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
507 }
508 }
509 }
510 }
511 // Catch clauses for all errors setting exc_type and exc_msg
512 OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
513
514 if (exc_type != ExceptionType::NONE) {
515 const std::string msg = "Compute initial well solution for new wells failed. Continue with zero initial rates";
516 local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
517 }
518
519 logAndCheckForExceptionsAndThrow(local_deferredLogger,
520 exc_type, "beginTimeStep() failed: " + exc_msg, terminal_output_, comm);
521
522 }
523
524 template<typename TypeTag>
525 void
526 BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx,
527 const double simulationTime,
528 DeferredLogger& deferred_logger)
529 {
530 for (const std::string& well_name : this->getWellsForTesting(timeStepIdx, simulationTime)) {
531 const Well& wellEcl = schedule().getWell(well_name, timeStepIdx);
532 if (wellEcl.getStatus() == Well::Status::SHUT)
533 continue;
534
535 WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
536 // some preparation before the well can be used
537 well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_, true);
538
539 double well_efficiency_factor = wellEcl.getEfficiencyFactor();
540 WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx),
541 schedule(), timeStepIdx, well_efficiency_factor);
542
543 well->setWellEfficiencyFactor(well_efficiency_factor);
544 well->setVFPProperties(vfp_properties_.get());
545 well->setGuideRate(&guideRate_);
546
547 well->wellTesting(ebosSimulator_, simulationTime, this->wellState(), this->groupState(), wellTestState(), deferred_logger);
548 }
549 }
550
551
552
553
554
555 // called at the end of a report step
556 template<typename TypeTag>
557 void
558 BlackoilWellModel<TypeTag>::
559 endReportStep()
560 {
561 // Clear the communication data structures for above values.
562 for (auto&& pinfo : this->local_parallel_well_info_)
563 {
564 pinfo.get().clear();
565 }
566 }
567
568
569
570
571
572 // called at the end of a report step
573 template<typename TypeTag>
574 const SimulatorReportSingle&
575 BlackoilWellModel<TypeTag>::
576 lastReport() const {return last_report_; }
577
578
579
580
581
582 // called at the end of a time step
583 template<typename TypeTag>
584 void
585 BlackoilWellModel<TypeTag>::
586 timeStepSucceeded(const double simulationTime, const double dt)
587 {
588 this->closed_this_step_.clear();
589
590 // time step is finished and we are not any more at the beginning of an report step
591 report_step_starts_ = false;
592 const int reportStepIdx = ebosSimulator_.episodeIndex();
593
594 DeferredLogger local_deferredLogger;
595 for (const auto& well : well_container_) {
596 if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
597 well->updateWaterThroughput(dt, this->wellState());
598 }
599 }
600
601 if (Indices::waterEnabled) {
602 this->updateFiltrationParticleVolume(dt, FluidSystem::waterPhaseIdx);
603 }
604
605 // at the end of the time step, updating the inj_multiplier saved in WellState for later use
606 this->updateInjMult(local_deferredLogger);
607
608 // report well switching
609 for (const auto& well : well_container_) {
610 well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);
611 }
612 // report group switching
613 if (terminal_output_) {
614
615 for (const auto& [name, to] : switched_prod_groups_) {
616 const Group::ProductionCMode& oldControl = this->prevWGState().group_state.production_control(name);
617 std::string from = Group::ProductionCMode2String(oldControl);
618 if (to != from) {
619 std::string msg = " Production Group " + name
620 + " control mode changed from ";
621 msg += from;
622 msg += " to " + to;
623 local_deferredLogger.info(msg);
624 }
625 }
626 for (const auto& [key, to] : switched_inj_groups_) {
627 const std::string& name = key.first;
628 const Opm::Phase& phase = key.second;
629
630 const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(name, phase);
631 std::string from = Group::InjectionCMode2String(oldControl);
632 if (to != from) {
633 std::string msg = " Injection Group " + name
634 + " control mode changed from ";
635 msg += from;
636 msg += " to " + to;
637 local_deferredLogger.info(msg);
638 }
639 }
640 }
641
642 // update the rate converter with current averages pressures etc in
643 rateConverter_->template defineState<ElementContext>(ebosSimulator_);
644
645 // calculate the well potentials
646 try {
647 updateWellPotentials(reportStepIdx,
648 /*onlyAfterEvent*/false,
649 ebosSimulator_.vanguard().summaryConfig(),
650 local_deferredLogger);
651 } catch ( std::runtime_error& e ) {
652 const std::string msg = "A zero well potential is returned for output purposes. ";
653 local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
654 }
655
656 updateWellTestState(simulationTime, wellTestState());
657
658 // check group sales limits at the end of the timestep
659 const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx);
660 checkGEconLimits(
661 fieldGroup, simulationTime, ebosSimulator_.episodeIndex(), local_deferredLogger);
662 checkGconsaleLimits(fieldGroup, this->wellState(),
663 ebosSimulator_.episodeIndex(), local_deferredLogger);
664
665 this->calculateProductivityIndexValues(local_deferredLogger);
666
667 this->commitWGState();
668
669 const Opm::Parallel::Communication& comm = grid().comm();
670 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
671 if (terminal_output_) {
672 global_deferredLogger.logMessages();
673 }
674
675 //reporting output temperatures
676 this->computeWellTemperature();
677 }
678
679
680 template<typename TypeTag>
681 void
682 BlackoilWellModel<TypeTag>::
683 computeTotalRatesForDof(RateVector& rate,
684 unsigned elemIdx) const
685 {
686 rate = 0;
687
688 if (!is_cell_perforated_[elemIdx])
689 return;
690
691 for (const auto& well : well_container_)
692 well->addCellRates(rate, elemIdx);
693 }
694
695
696 template<typename TypeTag>
697 template <class Context>
698 void
699 BlackoilWellModel<TypeTag>::
700 computeTotalRatesForDof(RateVector& rate,
701 const Context& context,
702 unsigned spaceIdx,
703 unsigned timeIdx) const
704 {
705 rate = 0;
706 int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
707
708 if (!is_cell_perforated_[elemIdx])
709 return;
710
711 for (const auto& well : well_container_)
712 well->addCellRates(rate, elemIdx);
713 }
714
715
716
717 template<typename TypeTag>
718 void
719 BlackoilWellModel<TypeTag>::
720 initializeWellState(const int timeStepIdx)
721 {
722 std::vector<double> cellPressures(this->local_num_cells_, 0.0);
723 ElementContext elemCtx(ebosSimulator_);
724
725 const auto& gridView = ebosSimulator_.vanguard().gridView();
726
727 OPM_BEGIN_PARALLEL_TRY_CATCH();
728 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
729 elemCtx.updatePrimaryStencil(elem);
730 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
731
732 const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState();
733 // copy of get perfpressure in Standard well except for value
734 double& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)];
735 if (Indices::oilEnabled) {
736 perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
737 } else if (Indices::waterEnabled) {
738 perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
739 } else {
740 perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
741 }
742 }
743 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ", ebosSimulator_.vanguard().grid().comm());
744
745 this->wellState().init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx,
746 &this->prevWellState(), well_perf_data_,
747 this->summaryState());
748 }
749
750
751
752
753
754 template<typename TypeTag>
755 void
756 BlackoilWellModel<TypeTag>::
757 createWellContainer(const int time_step)
758 {
759 DeferredLogger local_deferredLogger;
760
761 const int nw = numLocalWells();
762
763 well_container_.clear();
764
765 if (nw > 0) {
766 well_container_.reserve(nw);
767
768 for (int w = 0; w < nw; ++w) {
769 const Well& well_ecl = wells_ecl_[w];
770
771 if (!well_ecl.hasConnections()) {
772 // No connections in this well. Nothing to do.
773 continue;
774 }
775
776 const std::string& well_name = well_ecl.name();
777 const auto well_status = this->schedule()
778 .getWell(well_name, time_step).getStatus();
779
780 if ((well_ecl.getStatus() == Well::Status::SHUT) ||
781 (well_status == Well::Status::SHUT))
782 {
783 // Due to ACTIONX the well might have been closed behind our back.
784 if (well_ecl.getStatus() != Well::Status::SHUT) {
785 this->closed_this_step_.insert(well_name);
786 this->wellState().shutWell(w);
787 }
788
789 continue;
790 }
791
792 // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
793 if (this->wellTestState().well_is_closed(well_name)) {
794 // TODO: more checking here, to make sure this standard more specific and complete
795 // maybe there is some WCON keywords will not open the well
796 auto& events = this->wellState().well(w).events;
797 if (events.hasEvent(WellState::event_mask)) {
798 if (wellTestState().lastTestTime(well_name) == ebosSimulator_.time()) {
799 // The well was shut this timestep, we are most likely retrying
800 // a timestep without the well in question, after it caused
801 // repeated timestep cuts. It should therefore not be opened,
802 // even if it was new or received new targets this report step.
803 events.clearEvent(WellState::event_mask);
804 } else {
805 wellTestState().open_well(well_name);
806 wellTestState().open_completions(well_name);
807 }
808 }
809 }
810
811 // TODO: should we do this for all kinds of closing reasons?
812 // something like wellTestState().hasWell(well_name)?
813 bool wellIsStopped = false;
814 if (wellTestState().well_is_closed(well_name))
815 {
816 if (well_ecl.getAutomaticShutIn()) {
817 // shut wells are not added to the well container
818 this->wellState().shutWell(w);
819 continue;
820 } else {
821 if (!well_ecl.getAllowCrossFlow()) {
822 // stopped wells where cross flow is not allowed
823 // are not added to the well container
824 this->wellState().shutWell(w);
825 continue;
826 }
827 // stopped wells are added to the container but marked as stopped
828 this->wellState().stopWell(w);
829 wellIsStopped = true;
830 }
831 }
832
833 // If a production well disallows crossflow and its
834 // (prediction type) rate control is zero, then it is effectively shut.
835 if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) {
836 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
837 const auto prod_controls = well_ecl.productionControls(summaryState);
838
839 auto is_zero = [](const double x)
840 {
841 return std::isfinite(x) && !std::isnormal(x);
842 };
843
844 bool zero_rate_control = false;
845 switch (prod_controls.cmode) {
846 case Well::ProducerCMode::ORAT:
847 zero_rate_control = is_zero(prod_controls.oil_rate);
848 break;
849
850 case Well::ProducerCMode::WRAT:
851 zero_rate_control = is_zero(prod_controls.water_rate);
852 break;
853
854 case Well::ProducerCMode::GRAT:
855 zero_rate_control = is_zero(prod_controls.gas_rate);
856 break;
857
858 case Well::ProducerCMode::LRAT:
859 zero_rate_control = is_zero(prod_controls.liquid_rate);
860 break;
861
862 case Well::ProducerCMode::RESV:
863 zero_rate_control = is_zero(prod_controls.resv_rate);
864 break;
865 default:
866 // Might still have zero rate controls, but is pressure controlled.
867 zero_rate_control = false;
868 break;
869 }
870
871 if (zero_rate_control) {
872 // Treat as shut, do not add to container.
873 local_deferredLogger.info(" Well shut due to zero rate control and disallowing crossflow: " + well_ecl.name());
874 this->wellState().shutWell(w);
875 continue;
876 }
877 }
878
879 if (well_status == Well::Status::STOP) {
880 this->wellState().stopWell(w);
881 wellIsStopped = true;
882 }
883
884 well_container_.emplace_back(this->createWellPointer(w, time_step));
885
886 if (wellIsStopped)
887 well_container_.back()->stopWell();
888 }
889 }
890
891 // Collect log messages and print.
892
893 const Opm::Parallel::Communication& comm = grid().comm();
894 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
895 if (terminal_output_) {
896 global_deferredLogger.logMessages();
897 }
898
899 well_container_generic_.clear();
900 for (auto& w : well_container_)
901 well_container_generic_.push_back(w.get());
902
903 const auto& network = schedule()[time_step].network();
904 if (network.active() && !this->node_pressures_.empty()) {
905 for (auto& well: well_container_generic_) {
906 // Producers only, since we so far only support the
907 // "extended" network model (properties defined by
908 // BRANPROP and NODEPROP) which only applies to producers.
909 if (well->isProducer()) {
910 const auto it = node_pressures_.find(well->wellEcl().groupName());
911 if (it != node_pressures_.end()) {
912 // The well belongs to a group which has a network nodal pressure,
913 // set the dynamic THP constraint based on the network nodal pressure
914 const double nodal_pressure = it->second;
915 well->setDynamicThpLimit(nodal_pressure);
916 }
917 }
918 }
919 }
920
921 this->registerOpenWellsForWBPCalculation();
922 }
923
924
925
926
927
928 template <typename TypeTag>
929 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
930 BlackoilWellModel<TypeTag>::
931 createWellPointer(const int wellID, const int time_step) const
932 {
933 const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
934
935 if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
936 return this->template createTypedWellPointer<StandardWell<TypeTag>>(wellID, time_step);
937 }
938 else {
939 return this->template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, time_step);
940 }
941 }
942
943
944
945
946
947 template <typename TypeTag>
948 template <typename WellType>
949 std::unique_ptr<WellType>
950 BlackoilWellModel<TypeTag>::
951 createTypedWellPointer(const int wellID, const int time_step) const
952 {
953 // Use the pvtRegionIdx from the top cell
954 const auto& perf_data = this->well_perf_data_[wellID];
955
956 // Cater for case where local part might have no perforations.
957 const auto pvtreg = perf_data.empty()
958 ? 0 : pvt_region_idx_[perf_data.front().cell_index];
959
960 const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get();
961 const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
962
963 return std::make_unique<WellType>(this->wells_ecl_[wellID],
964 parallel_well_info,
965 time_step,
966 this->param_,
967 *this->rateConverter_,
968 global_pvtreg,
969 this->numComponents(),
970 this->numPhases(),
971 wellID,
972 perf_data);
973 }
974
975
976
977
978
979 template<typename TypeTag>
980 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
981 BlackoilWellModel<TypeTag>::
982 createWellForWellTest(const std::string& well_name,
983 const int report_step,
984 DeferredLogger& deferred_logger) const
985 {
986 // Finding the location of the well in wells_ecl
987 const int nw_wells_ecl = wells_ecl_.size();
988 int index_well_ecl = 0;
989 for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
990 if (well_name == wells_ecl_[index_well_ecl].name()) {
991 break;
992 }
993 }
994 // It should be able to find in wells_ecl.
995 if (index_well_ecl == nw_wells_ecl) {
996 OPM_DEFLOG_THROW(std::logic_error,
997 fmt::format("Could not find well {} in wells_ecl ", well_name),
998 deferred_logger);
999 }
1000
1001 return this->createWellPointer(index_well_ecl, report_step);
1002 }
1003
1004
1005
1006 template<typename TypeTag>
1007 void
1008 BlackoilWellModel<TypeTag>::
1009 balanceNetwork(DeferredLogger& deferred_logger) {
1010 const double dt = this->ebosSimulator_.timeStepSize();
1011 // TODO: should we also have the group and network backed-up here in case the solution did not get converged?
1012 auto& well_state = this->wellState();
1013 constexpr std::size_t max_iter = 100;
1014 bool converged = false;
1015 std::size_t iter = 0;
1016 bool changed_well_group = false;
1017 do {
1018 changed_well_group = updateWellControlsAndNetwork(true, dt, deferred_logger);
1019 assembleWellEqWithoutIteration(dt, deferred_logger);
1020 converged = this->getWellConvergence(this->B_avg_, true).converged() && !changed_well_group;
1021 if (converged) {
1022 break;
1023 }
1024 ++iter;
1025 for (auto& well : this->well_container_) {
1026 const auto& summary_state = this->ebosSimulator_.vanguard().summaryState();
1027 well->solveEqAndUpdateWellState(summary_state, well_state, deferred_logger);
1028 }
1029 this->initPrimaryVariablesEvaluation();
1030 } while (iter < max_iter);
1031
1032 if (!converged) {
1033 const std::string msg = fmt::format("balanceNetwork did not get converged with {} iterations, and unconverged "
1034 "network balance result will be used", max_iter);
1035 deferred_logger.warning(msg);
1036 } else {
1037 const std::string msg = fmt::format("balanceNetwork get converged with {} iterations", iter);
1038 deferred_logger.debug(msg);
1039 }
1040 }
1041
1042
1043
1044
1045 template<typename TypeTag>
1046 void
1047 BlackoilWellModel<TypeTag>::
1048 assemble(const int iterationIdx,
1049 const double dt)
1050 {
1051
1052 DeferredLogger local_deferredLogger;
1053 if (this->glift_debug) {
1054 const std::string msg = fmt::format(
1055 "assemble() : iteration {}" , iterationIdx);
1056 gliftDebug(msg, local_deferredLogger);
1057 }
1058 last_report_ = SimulatorReportSingle();
1059 Dune::Timer perfTimer;
1060 perfTimer.start();
1061
1062 {
1063 const int episodeIdx = ebosSimulator_.episodeIndex();
1064 const auto& network = schedule()[episodeIdx].network();
1065 if ( !wellsActive() && !network.active() ) {
1066 return;
1067 }
1068 }
1069
1070 if (iterationIdx == 0 && wellsActive()) {
1071 // try-catch is needed here as updateWellControls
1072 // contains global communication and has either to
1073 // be reached by all processes or all need to abort
1074 // before.
1075 OPM_BEGIN_PARALLEL_TRY_CATCH();
1076 {
1077 calculateExplicitQuantities(local_deferredLogger);
1078 prepareTimeStep(local_deferredLogger);
1079 }
1080 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1081 "assemble() failed (It=0): ",
1082 terminal_output_, grid().comm());
1083 }
1084
1085 const bool well_group_control_changed = updateWellControlsAndNetwork(false, dt, local_deferredLogger);
1086
1087 // even when there is no wells active, the network nodal pressure still need to be updated through updateWellControlsAndNetwork()
1088 // but there is no need to assemble the well equations
1089 if ( ! wellsActive() ) {
1090 return;
1091 }
1092
1093 assembleWellEqWithoutIteration(dt, local_deferredLogger);
1094
1095 // if group or well control changes we don't consider the
1096 // case converged
1097 last_report_.well_group_control_changed = well_group_control_changed;
1098 last_report_.assemble_time_well += perfTimer.stop();
1099 }
1100
1101
1102
1103
1104 template<typename TypeTag>
1105 bool
1106 BlackoilWellModel<TypeTag>::
1107 updateWellControlsAndNetwork(const bool mandatory_network_balance, const double dt, DeferredLogger& local_deferredLogger)
1108 {
1109 // not necessarily that we always need to update once of the network solutions
1110 bool do_network_update = true;
1111 bool well_group_control_changed = false;
1112 // after certain number of the iterations, we use relaxed tolerance for the network update
1113 const std::size_t iteration_to_relax = param_.network_max_strict_iterations_;
1114 // after certain number of the iterations, we terminate
1115 const std::size_t max_iteration = param_.network_max_iterations_;
1116 std::size_t network_update_iteration = 0;
1117 while (do_network_update) {
1118 if (terminal_output_ && (network_update_iteration == iteration_to_relax) ) {
1119 local_deferredLogger.info(" we begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations ");
1120 }
1121 const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
1122 std::tie(do_network_update, well_group_control_changed) =
1123 updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger);
1124 ++network_update_iteration;
1125
1126 if (network_update_iteration >= max_iteration ) {
1127 if (terminal_output_) {
1128 local_deferredLogger.info("maximum of " + std::to_string(max_iteration) + " iterations has been used, we stop the network update now. "
1129 "The simulation will continue with unconverged network results");
1130 }
1131 break;
1132 }
1133 }
1134 return well_group_control_changed;
1135 }
1136
1137
1138
1139
1140 template<typename TypeTag>
1141 std::pair<bool, bool>
1142 BlackoilWellModel<TypeTag>::
1143 updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
1144 const bool relax_network_tolerance,
1145 const double dt,
1146 DeferredLogger& local_deferredLogger)
1147 {
1148 auto [well_group_control_changed, more_network_update] =
1149 updateWellControls(mandatory_network_balance, local_deferredLogger, relax_network_tolerance);
1150
1151 bool alq_updated = false;
1152 OPM_BEGIN_PARALLEL_TRY_CATCH();
1153 {
1154 // Set the well primary variables based on the value of well solutions
1155 initPrimaryVariablesEvaluation();
1156
1157 alq_updated = maybeDoGasLiftOptimize(local_deferredLogger);
1158
1159 prepareWellsBeforeAssembling(dt, local_deferredLogger);
1160 }
1161 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "updateWellControlsAndNetworkIteration() failed: ",
1162 terminal_output_, grid().comm());
1163
1164 //update guide rates
1165 const int reportStepIdx = ebosSimulator_.episodeIndex();
1166 if (alq_updated || BlackoilWellModelGuideRates(*this).
1167 guideRateUpdateIsNeeded(reportStepIdx)) {
1168 const double simulationTime = ebosSimulator_.time();
1169 const auto& comm = ebosSimulator_.vanguard().grid().comm();
1170 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
1171 std::vector<double> pot(numPhases(), 0.0);
1172 const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
1173 WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
1174 this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
1175 }
1176
1177
1178 return {more_network_update, well_group_control_changed};
1179 }
1180
1181
1182
1183 template<typename TypeTag>
1184 void
1185 BlackoilWellModel<TypeTag>::
1186 assembleDomain([[maybe_unused]] const int iterationIdx,
1187 const double dt,
1188 const Domain& domain)
1189 {
1190 last_report_ = SimulatorReportSingle();
1191 Dune::Timer perfTimer;
1192 perfTimer.start();
1193
1194 {
1195 const int episodeIdx = ebosSimulator_.episodeIndex();
1196 const auto& network = schedule()[episodeIdx].network();
1197 if ( !wellsActive() && !network.active() ) {
1198 return;
1199 }
1200 }
1201
1202 // We assume that calculateExplicitQuantities() and
1203 // prepareTimeStep() have been called already for the entire
1204 // well model, so we do not need to do it here (when
1205 // iterationIdx is 0).
1206
1207 DeferredLogger local_deferredLogger;
1208 updateWellControlsDomain(local_deferredLogger, domain);
1209 initPrimaryVariablesEvaluationDomain(domain);
1210 assembleWellEqDomain(dt, domain, local_deferredLogger);
1211
1212 // TODO: errors here must be caught higher up, as this method is not called in parallel.
1213 // We will log errors on rank 0, but not other ranks for now.
1214 if (terminal_output_) {
1215 local_deferredLogger.logMessages();
1216 }
1217
1218 last_report_.converged = true;
1219 last_report_.assemble_time_well += perfTimer.stop();
1220 }
1221
1222
1223 template<typename TypeTag>
1224 bool
1225 BlackoilWellModel<TypeTag>::
1226 maybeDoGasLiftOptimize(DeferredLogger& deferred_logger)
1227 {
1228 bool do_glift_optimization = false;
1229 int num_wells_changed = 0;
1230 const double simulation_time = ebosSimulator_.time();
1231 const double min_wait = ebosSimulator_.vanguard().schedule().glo(ebosSimulator_.episodeIndex()).min_wait();
1232 // We only optimize if a min_wait time has past.
1233 // If all_newton is true we still want to optimize several times pr timestep
1234 // i.e. we also optimize if check simulation_time == last_glift_opt_time_
1235 // that is when the last_glift_opt_time is already updated with the current time step
1236 if ( simulation_time == last_glift_opt_time_ || simulation_time >= (last_glift_opt_time_ + min_wait)) {
1237 do_glift_optimization = true;
1238 last_glift_opt_time_ = simulation_time;
1239 }
1240
1241 if (do_glift_optimization) {
1242 GLiftOptWells glift_wells;
1243 GLiftProdWells prod_wells;
1244 GLiftWellStateMap state_map;
1245 // NOTE: To make GasLiftGroupInfo (see below) independent of the TypeTag
1246 // associated with *this (i.e. BlackoilWellModel<TypeTag>) we observe
1247 // that GasLiftGroupInfo's only dependence on *this is that it needs to
1248 // access the eclipse Wells in the well container (the eclipse Wells
1249 // themselves are independent of the TypeTag).
1250 // Hence, we extract them from the well container such that we can pass
1251 // them to the GasLiftGroupInfo constructor.
1252 GLiftEclWells ecl_well_map;
1253 initGliftEclWellMap(ecl_well_map);
1254 GasLiftGroupInfo group_info {
1255 ecl_well_map,
1256 ebosSimulator_.vanguard().schedule(),
1257 ebosSimulator_.vanguard().summaryState(),
1258 ebosSimulator_.episodeIndex(),
1259 ebosSimulator_.model().newtonMethod().numIterations(),
1260 phase_usage_,
1261 deferred_logger,
1262 this->wellState(),
1263 this->groupState(),
1264 ebosSimulator_.vanguard().grid().comm(),
1265 this->glift_debug
1266 };
1267 group_info.initialize();
1268 gasLiftOptimizationStage1(
1269 deferred_logger, prod_wells, glift_wells, group_info, state_map);
1270 gasLiftOptimizationStage2(
1271 deferred_logger, prod_wells, glift_wells, group_info, state_map,
1272 ebosSimulator_.episodeIndex());
1273 if (this->glift_debug) gliftDebugShowALQ(deferred_logger);
1274 num_wells_changed = glift_wells.size();
1275 }
1276 num_wells_changed = this->comm_.sum(num_wells_changed);
1277 return num_wells_changed > 0;
1278 }
1279
1280 template<typename TypeTag>
1281 void
1282 BlackoilWellModel<TypeTag>::
1283 gasLiftOptimizationStage1(DeferredLogger& deferred_logger,
1284 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
1285 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map)
1286 {
1287 auto comm = ebosSimulator_.vanguard().grid().comm();
1288 int num_procs = comm.size();
1289 // NOTE: Gas lift optimization stage 1 seems to be difficult
1290 // to do in parallel since the wells are optimized on different
1291 // processes and each process needs to know the current ALQ allocated
1292 // to each group it is a memeber of in order to check group limits and avoid
1293 // allocating more ALQ than necessary. (Surplus ALQ is removed in
1294 // stage 2). In stage1, as each well is adding ALQ, the current group ALQ needs
1295 // to be communicated to the other processes. But there is no common
1296 // synchronization point that all process will reach in the
1297 // runOptimizeLoop_() in GasLiftSingleWell.cpp.
1298 //
1299 // TODO: Maybe a better solution could be invented by distributing
1300 // wells according to certain parent groups. Then updated group rates
1301 // might not have to be communicated to the other processors.
1302
1303 // Currently, the best option seems to be to run this part sequentially
1304 // (not in parallel).
1305 //
1306 // TODO: The simplest approach seems to be if a) one process could take
1307 // ownership of all the wells (the union of all the wells in the
1308 // well_container_ of each process) then this process could do the
1309 // optimization, while the other processes could wait for it to
1310 // finish (e.g. comm.barrier()), or alternatively, b) if all
1311 // processes could take ownership of all the wells. Then there
1312 // would be no need for synchronization here..
1313 //
1314 for (int i = 0; i< num_procs; i++) {
1315 int num_rates_to_sync = 0; // communication variable
1316 GLiftSyncGroups groups_to_sync;
1317 if (comm.rank() == i) {
1318 // Run stage1: Optimize single wells while also checking group limits
1319 for (const auto& well : well_container_) {
1320 // NOTE: Only the wells in "group_info" needs to be optimized
1321 if (group_info.hasWell(well->name())) {
1322 gasLiftOptimizationStage1SingleWell(
1323 well.get(), deferred_logger, prod_wells, glift_wells,
1324 group_info, state_map, groups_to_sync
1325 );
1326 }
1327 }
1328 num_rates_to_sync = groups_to_sync.size();
1329 }
1330 num_rates_to_sync = comm.sum(num_rates_to_sync);
1331 if (num_rates_to_sync > 0) {
1332 std::vector<int> group_indexes;
1333 group_indexes.reserve(num_rates_to_sync);
1334 std::vector<double> group_alq_rates;
1335 group_alq_rates.reserve(num_rates_to_sync);
1336 std::vector<double> group_oil_rates;
1337 group_oil_rates.reserve(num_rates_to_sync);
1338 std::vector<double> group_gas_rates;
1339 group_gas_rates.reserve(num_rates_to_sync);
1340 std::vector<double> group_water_rates;
1341 group_water_rates.reserve(num_rates_to_sync);
1342 if (comm.rank() == i) {
1343 for (auto idx : groups_to_sync) {
1344 auto [oil_rate, gas_rate, water_rate, alq] = group_info.getRates(idx);
1345 group_indexes.push_back(idx);
1346 group_oil_rates.push_back(oil_rate);
1347 group_gas_rates.push_back(gas_rate);
1348 group_water_rates.push_back(water_rate);
1349 group_alq_rates.push_back(alq);
1350 }
1351 } else {
1352 group_indexes.resize(num_rates_to_sync);
1353 group_oil_rates.resize(num_rates_to_sync);
1354 group_gas_rates.resize(num_rates_to_sync);
1355 group_water_rates.resize(num_rates_to_sync);
1356 group_alq_rates.resize(num_rates_to_sync);
1357 }
1358#if HAVE_MPI
1359 EclMpiSerializer ser(comm);
1360 ser.broadcast(i, group_indexes, group_oil_rates,
1361 group_gas_rates, group_water_rates, group_alq_rates);
1362#endif
1363 if (comm.rank() != i) {
1364 for (int j=0; j<num_rates_to_sync; j++) {
1365 group_info.updateRate(group_indexes[j],
1366 group_oil_rates[j], group_gas_rates[j], group_water_rates[j], group_alq_rates[j]);
1367 }
1368 }
1369 if (this->glift_debug) {
1370 int counter = 0;
1371 if (comm.rank() == i) {
1372 counter = this->wellState().gliftGetDebugCounter();
1373 }
1374 counter = comm.sum(counter);
1375 if (comm.rank() != i) {
1376 this->wellState().gliftSetDebugCounter(counter);
1377 }
1378 }
1379 }
1380 }
1381 }
1382
1383 // NOTE: this method cannot be const since it passes this->wellState()
1384 // (see below) to the GasLiftSingleWell constructor which accepts WellState
1385 // as a non-const reference..
1386 template<typename TypeTag>
1387 void
1388 BlackoilWellModel<TypeTag>::
1389 gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
1390 DeferredLogger& deferred_logger,
1391 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
1392 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
1393 GLiftSyncGroups& sync_groups)
1394 {
1395 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
1396 std::unique_ptr<GasLiftSingleWell> glift
1397 = std::make_unique<GasLiftSingleWell>(
1398 *well, ebosSimulator_, summary_state,
1399 deferred_logger, this->wellState(), this->groupState(),
1400 group_info, sync_groups, this->comm_, this->glift_debug);
1401 auto state = glift->runOptimize(
1402 ebosSimulator_.model().newtonMethod().numIterations());
1403 if (state) {
1404 state_map.insert({well->name(), std::move(state)});
1405 glift_wells.insert({well->name(), std::move(glift)});
1406 return;
1407 }
1408 prod_wells.insert({well->name(), well});
1409 }
1410
1411
1412 template<typename TypeTag>
1413 void
1414 BlackoilWellModel<TypeTag>::
1415 initGliftEclWellMap(GLiftEclWells &ecl_well_map)
1416 {
1417 for ( const auto& well: well_container_ ) {
1418 ecl_well_map.try_emplace(
1419 well->name(), &(well->wellEcl()), well->indexOfWell());
1420 }
1421 }
1422
1423
1424 template<typename TypeTag>
1425 void
1426 BlackoilWellModel<TypeTag>::
1427 assembleWellEq(const double dt, DeferredLogger& deferred_logger)
1428 {
1429 for (auto& well : well_container_) {
1430 well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1431 }
1432 }
1433
1434
1435 template<typename TypeTag>
1436 void
1437 BlackoilWellModel<TypeTag>::
1438 assembleWellEqDomain(const double dt, const Domain& domain, DeferredLogger& deferred_logger)
1439 {
1440 for (auto& well : well_container_) {
1441 if (well_domain_.at(well->name()) == domain.index) {
1442 well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1443 }
1444 }
1445 }
1446
1447
1448 template<typename TypeTag>
1449 void
1450 BlackoilWellModel<TypeTag>::
1451 prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger)
1452 {
1453 for (auto& well : well_container_) {
1454 well->prepareWellBeforeAssembling(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1455 }
1456 }
1457
1458
1459 template<typename TypeTag>
1460 void
1461 BlackoilWellModel<TypeTag>::
1462 assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger)
1463 {
1464 for (auto& well: well_container_) {
1465 well->assembleWellEqWithoutIteration(ebosSimulator_, dt, this->wellState(), this->groupState(),
1466 deferred_logger);
1467 }
1468 }
1469
1470
1471 template<typename TypeTag>
1472 void
1473 BlackoilWellModel<TypeTag>::
1474 apply(BVector& r) const
1475 {
1476 for (auto& well : well_container_) {
1477 well->apply(r);
1478 }
1479 }
1480
1481
1482 // Ax = A x - C D^-1 B x
1483 template<typename TypeTag>
1484 void
1485 BlackoilWellModel<TypeTag>::
1486 apply(const BVector& x, BVector& Ax) const
1487 {
1488 for (auto& well : well_container_) {
1489 well->apply(x, Ax);
1490 }
1491 }
1492
1493 template<typename TypeTag>
1494 void
1495 BlackoilWellModel<TypeTag>::
1496 getWellContributions(WellContributions& wellContribs) const
1497 {
1498 // prepare for StandardWells
1499 wellContribs.setBlockSize(StandardWell<TypeTag>::Indices::numEq, StandardWell<TypeTag>::numStaticWellEq);
1500
1501 for(unsigned int i = 0; i < well_container_.size(); i++){
1502 auto& well = well_container_[i];
1503 std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1504 if (derived) {
1505 wellContribs.addNumBlocks(derived->linSys().getNumBlocks());
1506 }
1507 }
1508
1509 // allocate memory for data from StandardWells
1510 wellContribs.alloc();
1511
1512 for(unsigned int i = 0; i < well_container_.size(); i++){
1513 auto& well = well_container_[i];
1514 // maybe WellInterface could implement addWellContribution()
1515 auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag>>(well);
1516 if (derived_std) {
1517 derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs);
1518 } else {
1519 auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
1520 if (derived_ms) {
1521 derived_ms->linSys().extract(wellContribs);
1522 } else {
1523 OpmLog::warning("Warning unknown type of well");
1524 }
1525 }
1526 }
1527 }
1528
1529 // Ax = Ax - alpha * C D^-1 B x
1530 template<typename TypeTag>
1531 void
1532 BlackoilWellModel<TypeTag>::
1533 applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
1534 {
1535 if (this->well_container_.empty()) {
1536 return;
1537 }
1538
1539 if( scaleAddRes_.size() != Ax.size() ) {
1540 scaleAddRes_.resize( Ax.size() );
1541 }
1542
1543 scaleAddRes_ = 0.0;
1544 // scaleAddRes_ = - C D^-1 B x
1545 apply( x, scaleAddRes_ );
1546 // Ax = Ax + alpha * scaleAddRes_
1547 Ax.axpy( alpha, scaleAddRes_ );
1548 }
1549
1550 template<typename TypeTag>
1551 void
1552 BlackoilWellModel<TypeTag>::
1553 addWellContributions(SparseMatrixAdapter& jacobian) const
1554 {
1555 for ( const auto& well: well_container_ ) {
1556 well->addWellContributions(jacobian);
1557 }
1558 }
1559
1560 template<typename TypeTag>
1561 void
1562 BlackoilWellModel<TypeTag>::
1563 addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const
1564 {
1565 int nw = this->numLocalWellsEnd();
1566 int rdofs = local_num_cells_;
1567 for ( int i = 0; i < nw; i++ ){
1568 int wdof = rdofs + i;
1569 jacobian[wdof][wdof] = 1.0;// better scaling ?
1570 }
1571
1572 for ( const auto& well : well_container_ ) {
1573 well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
1574 }
1575 }
1576
1577 template <typename TypeTag>
1578 void BlackoilWellModel<TypeTag>::
1579 addReservoirSourceTerms(GlobalEqVector& residual,
1580 std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress) const
1581 {
1582 // NB this loop may write multiple times to the same element
1583 // if a cell is perforated by more than one well, so it should
1584 // not be OpenMP-parallelized.
1585 for (const auto& well : well_container_) {
1586 if (!well->isOperableAndSolvable() && !well->wellIsStopped()) {
1587 continue;
1588 }
1589 const auto& cells = well->cells();
1590 const auto& rates = well->connectionRates();
1591 for (unsigned perfIdx = 0; perfIdx < rates.size(); ++perfIdx) {
1592 unsigned cellIdx = cells[perfIdx];
1593 auto rate = rates[perfIdx];
1594 rate *= -1.0;
1595 VectorBlockType res(0.0);
1596 using MatrixBlockType = typename SparseMatrixAdapter::MatrixBlock;
1597 MatrixBlockType bMat(0.0);
1598 ebosSimulator_.model().linearizer().setResAndJacobi(res, bMat, rate);
1599 residual[cellIdx] += res;
1600 *diagMatAddress[cellIdx] += bMat;
1601 }
1602 }
1603 }
1604
1605
1606 template<typename TypeTag>
1607 void
1608 BlackoilWellModel<TypeTag>::
1609 addWellPressureEquationsStruct(PressureMatrix& jacobian) const
1610 {
1611 int nw = this->numLocalWellsEnd();
1612 int rdofs = local_num_cells_;
1613 for(int i=0; i < nw; i++){
1614 int wdof = rdofs + i;
1615 jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
1616 }
1617 std::vector<std::vector<int>> wellconnections = getMaxWellConnections();
1618 for(int i=0; i < nw; i++){
1619 const auto& perfcells = wellconnections[i];
1620 for(int perfcell : perfcells){
1621 int wdof = rdofs + i;
1622 jacobian.entry(wdof,perfcell) = 0.0;
1623 jacobian.entry(perfcell, wdof) = 0.0;
1624 }
1625 }
1626 }
1627
1628
1629 template<typename TypeTag>
1630 void
1631 BlackoilWellModel<TypeTag>::
1632 recoverWellSolutionAndUpdateWellState(const BVector& x)
1633 {
1634 DeferredLogger local_deferredLogger;
1635 OPM_BEGIN_PARALLEL_TRY_CATCH();
1636 {
1637 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
1638 for (auto& well : well_container_) {
1639 well->recoverWellSolutionAndUpdateWellState(summary_state, x, this->wellState(), local_deferredLogger);
1640 }
1641 }
1642 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1643 "recoverWellSolutionAndUpdateWellState() failed: ",
1644 terminal_output_, ebosSimulator_.vanguard().grid().comm());
1645 }
1646
1647
1648 template<typename TypeTag>
1649 void
1650 BlackoilWellModel<TypeTag>::
1651 recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const Domain& domain)
1652 {
1653 // Note: no point in trying to do a parallel gathering
1654 // try/catch here, as this function is not called in
1655 // parallel but for each individual domain of each rank.
1656 DeferredLogger local_deferredLogger;
1657 const auto& summary_state = this->ebosSimulator_.vanguard().summaryState();
1658 for (auto& well : well_container_) {
1659 if (well_domain_.at(well->name()) == domain.index) {
1660 well->recoverWellSolutionAndUpdateWellState(summary_state, x,
1661 this->wellState(),
1662 local_deferredLogger);
1663 }
1664 }
1665 // TODO: avoid losing the logging information that could
1666 // be stored in the local_deferredlogger in a parallel case.
1667 if (terminal_output_) {
1668 local_deferredLogger.logMessages();
1669 }
1670 }
1671
1672
1673 template<typename TypeTag>
1674 void
1675 BlackoilWellModel<TypeTag>::
1676 initPrimaryVariablesEvaluation() const
1677 {
1678 for (auto& well : well_container_) {
1679 well->initPrimaryVariablesEvaluation();
1680 }
1681 }
1682
1683
1684 template<typename TypeTag>
1685 void
1686 BlackoilWellModel<TypeTag>::
1687 initPrimaryVariablesEvaluationDomain(const Domain& domain) const
1688 {
1689 for (auto& well : well_container_) {
1690 if (well_domain_.at(well->name()) == domain.index) {
1691 well->initPrimaryVariablesEvaluation();
1692 }
1693 }
1694 }
1695
1696
1697
1698
1699
1700
1701 template<typename TypeTag>
1702 ConvergenceReport
1703 BlackoilWellModel<TypeTag>::
1704 getDomainWellConvergence(const Domain& domain,
1705 const std::vector<Scalar>& B_avg) const
1706 {
1707 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
1708 const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1709 const bool relax_tolerance = iterationIdx > param_.strict_outer_iter_wells_;
1710
1711 Opm::DeferredLogger local_deferredLogger;
1712 ConvergenceReport local_report;
1713 for (const auto& well : well_container_) {
1714 if ((well_domain_.at(well->name()) == domain.index)) {
1715 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
1716 local_report += well->getWellConvergence(summary_state,
1717 this->wellState(),
1718 B_avg,
1719 local_deferredLogger,
1720 relax_tolerance);
1721 } else {
1722 ConvergenceReport report;
1723 using CR = ConvergenceReport;
1724 report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
1725 local_report += report;
1726 }
1727 }
1728 }
1729
1730 // This function will be called once for each domain on a process,
1731 // and the number of domains on a process will vary, so there is
1732 // no way to communicate here. There is also no need, as a domain
1733 // is local to a single process in our current approach.
1734 // Therefore there is no call to gatherDeferredLogger() or to
1735 // gatherConvergenceReport() below. However, as of now, log messages
1736 // on non-output ranks will be lost here.
1737 // TODO: create the DeferredLogger outside this scope to enable it
1738 // to be gathered and printed on the output rank.
1739 Opm::DeferredLogger global_deferredLogger = local_deferredLogger;
1740 ConvergenceReport report = local_report;
1741 if (terminal_output_) {
1742 global_deferredLogger.logMessages();
1743 }
1744
1745 // Log debug messages for NaN or too large residuals.
1746 // TODO: This will as of now only be logged on the output rank.
1747 // In the similar code in getWellConvergence(), all ranks will be
1748 // at the same spot, that does not hold for this per-domain function.
1749 if (terminal_output_) {
1750 for (const auto& f : report.wellFailures()) {
1751 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
1752 OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1753 } else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
1754 OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1755 }
1756 }
1757 }
1758 return report;
1759 }
1760
1761
1762
1763
1764
1765 template<typename TypeTag>
1766 ConvergenceReport
1767 BlackoilWellModel<TypeTag>::
1768 getWellConvergence(const std::vector<Scalar>& B_avg, bool checkWellGroupControls) const
1769 {
1770
1771 DeferredLogger local_deferredLogger;
1772 // Get global (from all processes) convergence report.
1773 ConvergenceReport local_report;
1774 const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1775 for (const auto& well : well_container_) {
1776 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
1777 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
1778 local_report += well->getWellConvergence(
1779 summary_state, this->wellState(), B_avg, local_deferredLogger,
1780 iterationIdx > param_.strict_outer_iter_wells_);
1781 } else {
1782 ConvergenceReport report;
1783 using CR = ConvergenceReport;
1784 report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
1785 local_report += report;
1786 }
1787 }
1788
1789 const Opm::Parallel::Communication comm = grid().comm();
1790 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
1791 ConvergenceReport report = gatherConvergenceReport(local_report, comm);
1792
1793 // the well_group_control_changed info is already communicated
1794 if (checkWellGroupControls) {
1795 report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed);
1796 }
1797
1798 if (terminal_output_) {
1799 global_deferredLogger.logMessages();
1800 }
1801 // Log debug messages for NaN or too large residuals.
1802 if (terminal_output_) {
1803 for (const auto& f : report.wellFailures()) {
1804 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
1805 OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1806 } else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
1807 OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1808 }
1809 }
1810 }
1811 return report;
1812 }
1813
1814
1815
1816
1817
1818 template<typename TypeTag>
1819 void
1822 {
1823 // TODO: checking isOperableAndSolvable() ?
1824 for (auto& well : well_container_) {
1825 well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), deferred_logger);
1826 }
1827 }
1828
1829
1830
1831
1832
1833 template<typename TypeTag>
1834 std::pair<bool, bool>
1837 {
1838 const int episodeIdx = ebosSimulator_.episodeIndex();
1839 const auto& network = schedule()[episodeIdx].network();
1840 if (!wellsActive() && !network.active()) {
1841 return {false, false};
1842 }
1843
1844 const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1845 const auto& comm = ebosSimulator_.vanguard().grid().comm();
1846 updateAndCommunicateGroupData(episodeIdx, iterationIdx);
1847
1848 // network related
1849 bool more_network_update = false;
1850 if (shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
1851 const auto local_network_imbalance = updateNetworkPressures(episodeIdx);
1852 const double network_imbalance = comm.max(local_network_imbalance);
1853 const auto& balance = schedule()[episodeIdx].network_balance();
1854 constexpr double relaxtion_factor = 10.0;
1855 const double tolerance = relax_network_tolerance ? relaxtion_factor * balance.pressure_tolerance() : balance.pressure_tolerance();
1856 more_network_update = network_imbalance > tolerance;
1857 }
1858
1859 bool changed_well_group = false;
1860 // Check group individual constraints.
1861 const int nupcol = schedule()[episodeIdx].nupcol();
1862 // don't switch group control when iterationIdx > nupcol
1863 // to avoid oscilations between group controls
1864 if (iterationIdx <= nupcol) {
1865 const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
1866 changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
1867 }
1868 // Check wells' group constraints and communicate.
1869 bool changed_well_to_group = false;
1870 for (const auto& well : well_container_) {
1871 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
1872 const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1873 if (changed_well) {
1874 changed_well_to_group = changed_well || changed_well_to_group;
1875 }
1876 }
1877
1878 changed_well_to_group = comm.sum(changed_well_to_group);
1879 if (changed_well_to_group) {
1880 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1881 changed_well_group = true;
1882 }
1883
1884 // Check individual well constraints and communicate.
1885 bool changed_well_individual = false;
1886 for (const auto& well : well_container_) {
1887 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
1888 const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1889 if (changed_well) {
1890 changed_well_individual = changed_well || changed_well_individual;
1891 }
1892 }
1893 changed_well_individual = comm.sum(changed_well_individual);
1894 if (changed_well_individual) {
1895 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1896 changed_well_group = true;
1897 }
1898
1899 // update wsolvent fraction for REIN wells
1900 const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
1901 updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState());
1902
1903 return { changed_well_group, more_network_update };
1904 }
1905
1906 template<typename TypeTag>
1907 void
1908 BlackoilWellModel<TypeTag>::
1909 updateWellControlsDomain(DeferredLogger& deferred_logger, const Domain& domain)
1910 {
1911 if ( !wellsActive() ) return ;
1912
1913 // TODO: decide on and implement an approach to handling of
1914 // group controls, network and similar for domain solves.
1915
1916 // Check only individual well constraints and communicate.
1917 for (const auto& well : well_container_) {
1918 if (well_domain_.at(well->name()) == domain.index) {
1919 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
1920 well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1921 }
1922 }
1923 }
1924
1925
1926
1927
1928
1929 template <typename TypeTag>
1930 void
1931 BlackoilWellModel<TypeTag>::
1932 initializeWBPCalculationService()
1933 {
1934 this->wbpCalcMap_.clear();
1935 this->wbpCalcMap_.resize(this->wells_ecl_.size());
1936
1937 this->registerOpenWellsForWBPCalculation();
1938
1939 auto wellID = std::size_t{0};
1940 for (const auto& well : this->wells_ecl_) {
1941 this->wbpCalcMap_[wellID].wbpCalcIdx_ = this->wbpCalculationService_
1942 .createCalculator(well,
1943 this->local_parallel_well_info_[wellID],
1944 this->conn_idx_map_[wellID].local(),
1945 this->makeWellSourceEvaluatorFactory(wellID));
1946
1947 ++wellID;
1948 }
1949
1950 this->wbpCalculationService_.defineCommunication();
1951 }
1952
1953
1954
1955
1956
1957 template <typename TypeTag>
1958 data::WellBlockAveragePressures
1959 BlackoilWellModel<TypeTag>::
1960 computeWellBlockAveragePressures() const
1961 {
1962 auto wbpResult = data::WellBlockAveragePressures{};
1963
1964 using Calculated = PAvgCalculator::Result::WBPMode;
1965 using Output = data::WellBlockAvgPress::Quantity;
1966
1967 this->wbpCalculationService_.collectDynamicValues();
1968
1969 const auto numWells = this->wells_ecl_.size();
1970 for (auto wellID = 0*numWells; wellID < numWells; ++wellID) {
1971 const auto calcIdx = this->wbpCalcMap_[wellID].wbpCalcIdx_;
1972 const auto& well = this->wells_ecl_[wellID];
1973
1974 if (! well.hasRefDepth()) {
1975 // Can't perform depth correction without at least a
1976 // fall-back datum depth.
1977 continue;
1978 }
1979
1980 this->wbpCalculationService_
1981 .inferBlockAveragePressures(calcIdx, well.pavg(),
1982 this->gravity_,
1983 well.getWPaveRefDepth());
1984
1985 const auto& result = this->wbpCalculationService_
1986 .averagePressures(calcIdx);
1987
1988 auto& reported = wbpResult.values[well.name()];
1989
1990 reported[Output::WBP] = result.value(Calculated::WBP);
1991 reported[Output::WBP4] = result.value(Calculated::WBP4);
1992 reported[Output::WBP5] = result.value(Calculated::WBP5);
1993 reported[Output::WBP9] = result.value(Calculated::WBP9);
1994 }
1995
1996 return wbpResult;
1997 }
1998
1999
2000
2001
2002
2003 template <typename TypeTag>
2004 ParallelWBPCalculation::EvaluatorFactory
2005 BlackoilWellModel<TypeTag>::
2006 makeWellSourceEvaluatorFactory(const std::vector<Well>::size_type wellIdx) const
2007 {
2008 using Span = PAvgDynamicSourceData::SourceDataSpan<double>;
2009 using Item = typename Span::Item;
2010
2011 return [wellIdx, this]() -> ParallelWBPCalculation::Evaluator
2012 {
2013 if (! this->wbpCalcMap_[wellIdx].openWellIdx_.has_value()) {
2014 // Well is stopped/shut. Return evaluator for stopped wells.
2015 return []([[maybe_unused]] const int connIdx, Span sourceTerm)
2016 {
2017 // Well/connection is stopped/shut. Set all items to
2018 // zero.
2019
2020 sourceTerm
2021 .set(Item::Pressure , 0.0)
2022 .set(Item::PoreVol , 0.0)
2023 .set(Item::MixtureDensity, 0.0);
2024 };
2025 }
2026
2027 // Well is open. Return an evaluator for open wells/open connections.
2028 return [this, wellPtr = this->well_container_[*this->wbpCalcMap_[wellIdx].openWellIdx_].get()]
2029 (const int connIdx, Span sourceTerm)
2030 {
2031 // Note: The only item which actually matters for the WBP
2032 // calculation at the well reservoir connection level is the
2033 // mixture density. Set other items to zero.
2034
2035 const auto& connIdxMap =
2036 this->conn_idx_map_[wellPtr->indexOfWell()];
2037
2038 const auto rho = wellPtr->
2039 connectionDensity(connIdxMap.global(connIdx),
2040 connIdxMap.open(connIdx));
2041
2042 sourceTerm
2043 .set(Item::Pressure , 0.0)
2044 .set(Item::PoreVol , 0.0)
2045 .set(Item::MixtureDensity, rho);
2046 };
2047 };
2048 }
2049
2050
2051
2052
2053
2054 template <typename TypeTag>
2055 void
2056 BlackoilWellModel<TypeTag>::
2057 registerOpenWellsForWBPCalculation()
2058 {
2059 assert (this->wbpCalcMap_.size() == this->wells_ecl_.size());
2060
2061 for (auto& wbpCalc : this->wbpCalcMap_) {
2062 wbpCalc.openWellIdx_.reset();
2063 }
2064
2065 auto openWellIdx = typename std::vector<WellInterfacePtr>::size_type{0};
2066 for (const auto* openWell : this->well_container_generic_) {
2067 this->wbpCalcMap_[openWell->indexOfWell()].openWellIdx_ = openWellIdx++;
2068 }
2069 }
2070
2071
2072
2073
2074
2075 template<typename TypeTag>
2076 void
2077 BlackoilWellModel<TypeTag>::
2078 updateAndCommunicate(const int reportStepIdx,
2079 const int iterationIdx,
2080 DeferredLogger& deferred_logger)
2081 {
2082 updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
2083 // if a well or group change control it affects all wells that are under the same group
2084 for (const auto& well : well_container_) {
2085 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
2086 }
2087 updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
2088 }
2089
2090 template<typename TypeTag>
2091 bool
2092 BlackoilWellModel<TypeTag>::
2093 updateGroupControls(const Group& group,
2094 DeferredLogger& deferred_logger,
2095 const int reportStepIdx,
2096 const int iterationIdx)
2097 {
2098 bool changed = false;
2099 bool changed_hc = checkGroupHigherConstraints( group, deferred_logger, reportStepIdx);
2100 if (changed_hc) {
2101 changed = true;
2102 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
2103 }
2104 bool changed_individual =
2105 BlackoilWellModelConstraints(*this).
2106 updateGroupIndividualControl(group,
2107 reportStepIdx,
2108 this->switched_inj_groups_,
2109 this->switched_prod_groups_,
2110 this->groupState(),
2111 this->wellState(),
2112 deferred_logger);
2113 if (changed_individual) {
2114 changed = true;
2115 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
2116 }
2117 // call recursively down the group hierarchy
2118 for (const std::string& groupName : group.groups()) {
2119 bool changed_this = updateGroupControls( schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx,iterationIdx);
2120 changed = changed || changed_this;
2121 }
2122 return changed;
2123 }
2124
2125 template<typename TypeTag>
2126 void
2128 updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
2129 {
2131 for (const auto& well : well_container_) {
2132 const auto& wname = well->name();
2133 const auto wasClosed = wellTestState.well_is_closed(wname);
2134 well->checkWellOperability(ebosSimulator_, this->wellState(), local_deferredLogger);
2135 well->updateWellTestState(this->wellState().well(wname), simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
2136
2137 if (!wasClosed && wellTestState.well_is_closed(wname)) {
2138 this->closed_this_step_.insert(wname);
2139 }
2140 }
2141
2142 const Opm::Parallel::Communication comm = grid().comm();
2144 if (terminal_output_) {
2145 global_deferredLogger.logMessages();
2146 }
2147 }
2148
2149
2150 template<typename TypeTag>
2151 void
2154 std::string& exc_msg,
2155 ExceptionType::ExcEnum& exc_type,
2157 {
2158 const int np = numPhases();
2159 std::vector<double> potentials;
2160 const auto& well= well_container_[widx];
2161 try {
2162 well->computeWellPotentials(ebosSimulator_, well_state_copy, potentials, deferred_logger);
2163 }
2164 // catch all possible exception and store type and message.
2165 OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
2166 // Store it in the well state
2167 // potentials is resized and set to zero in the beginning of well->ComputeWellPotentials
2168 // and updated only if sucessfull. i.e. the potentials are zero for exceptions
2169 auto& ws = this->wellState().well(well->indexOfWell());
2170 for (int p = 0; p < np; ++p) {
2171 // make sure the potentials are positive
2172 ws.well_potentials[p] = std::max(0.0, potentials[p]);
2173 }
2174 }
2175
2176
2177
2178 template <typename TypeTag>
2179 void
2180 BlackoilWellModel<TypeTag>::
2181 calculateProductivityIndexValues(DeferredLogger& deferred_logger)
2182 {
2183 for (const auto& wellPtr : this->well_container_) {
2184 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2185 }
2186 }
2187
2188
2189
2190
2191
2192 template <typename TypeTag>
2193 void
2194 BlackoilWellModel<TypeTag>::
2195 calculateProductivityIndexValuesShutWells(const int reportStepIdx,
2196 DeferredLogger& deferred_logger)
2197 {
2198 // For the purpose of computing PI/II values, it is sufficient to
2199 // construct StandardWell instances only. We don't need to form
2200 // well objects that honour the 'isMultisegment()' flag of the
2201 // corresponding "this->wells_ecl_[shutWell]".
2202
2203 for (const auto& shutWell : this->local_shut_wells_) {
2204 if (!this->wells_ecl_[shutWell].hasConnections()) {
2205 // No connections in this well. Nothing to do.
2206 continue;
2207 }
2208
2209 auto wellPtr = this->template createTypedWellPointer
2210 <StandardWell<TypeTag>>(shutWell, reportStepIdx);
2211
2212 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
2213 this->local_num_cells_, this->B_avg_, true);
2214
2215 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2216 }
2217 }
2218
2219
2220
2221
2222
2223 template <typename TypeTag>
2224 void
2225 BlackoilWellModel<TypeTag>::
2226 calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
2227 DeferredLogger& deferred_logger)
2228 {
2229 wellPtr->updateProductivityIndex(this->ebosSimulator_,
2230 this->prod_index_calc_[wellPtr->indexOfWell()],
2231 this->wellState(),
2232 deferred_logger);
2233 }
2234
2235
2236
2237 template<typename TypeTag>
2238 void
2239 BlackoilWellModel<TypeTag>::
2240 prepareTimeStep(DeferredLogger& deferred_logger)
2241 {
2242 const bool network_rebalance_necessary = this->needRebalanceNetwork(ebosSimulator_.episodeIndex());
2243
2244 for (const auto& well : well_container_) {
2245 auto& events = this->wellState().well(well->indexOfWell()).events;
2246 if (events.hasEvent(WellState::event_mask)) {
2247 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
2248 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
2249 well->updatePrimaryVariables(summary_state, this->wellState(), deferred_logger);
2250 well->initPrimaryVariablesEvaluation();
2251 // There is no new well control change input within a report step,
2252 // so next time step, the well does not consider to have effective events anymore.
2253 events.clearEvent(WellState::event_mask);
2254 }
2255 // solve the well equation initially to improve the initial solution of the well model
2256 if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
2257 try {
2258 well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), deferred_logger);
2259 } catch (const std::exception& e) {
2260 const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the privious rates";
2261 deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
2262 }
2263 }
2264 // If we're using local well solves that include control switches, they also update
2265 // operability, so reset before main iterations begin
2266 well->resetWellOperability();
2267 }
2268 updatePrimaryVariables(deferred_logger);
2269
2270 if (network_rebalance_necessary) {
2271 // this is to obtain good network solution
2272 balanceNetwork(deferred_logger);
2273 }
2274 }
2275
2276 template<typename TypeTag>
2277 void
2278 BlackoilWellModel<TypeTag>::
2279 updateAverageFormationFactor()
2280 {
2281 std::vector< Scalar > B_avg(numComponents(), Scalar() );
2282 const auto& grid = ebosSimulator_.vanguard().grid();
2283 const auto& gridView = grid.leafGridView();
2284 ElementContext elemCtx(ebosSimulator_);
2285
2286 OPM_BEGIN_PARALLEL_TRY_CATCH();
2287 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2288 elemCtx.updatePrimaryStencil(elem);
2289 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
2290
2291 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
2292 const auto& fs = intQuants.fluidState();
2293
2294 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
2295 {
2296 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2297 continue;
2298 }
2299
2300 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2301 auto& B = B_avg[ compIdx ];
2302
2303 B += 1 / fs.invB(phaseIdx).value();
2304 }
2305 if constexpr (has_solvent_) {
2306 auto& B = B_avg[solventSaturationIdx];
2307 B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
2308 }
2309 }
2310 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updateAverageFormationFactor() failed: ", grid.comm())
2311
2312 // compute global average
2313 grid.comm().sum(B_avg.data(), B_avg.size());
2314 for (auto& bval : B_avg)
2315 {
2316 bval /= global_num_cells_;
2317 }
2318 B_avg_ = B_avg;
2319 }
2320
2321
2322
2323
2324
2325 template<typename TypeTag>
2326 void
2327 BlackoilWellModel<TypeTag>::
2328 updatePrimaryVariables(DeferredLogger& deferred_logger)
2329 {
2330 for (const auto& well : well_container_) {
2331 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
2332 well->updatePrimaryVariables(summary_state, this->wellState(), deferred_logger);
2333 }
2334 }
2335
2336 template<typename TypeTag>
2337 void
2338 BlackoilWellModel<TypeTag>::extractLegacyCellPvtRegionIndex_()
2339 {
2340 const auto& grid = ebosSimulator_.vanguard().grid();
2341 const auto& eclProblem = ebosSimulator_.problem();
2342 const unsigned numCells = grid.size(/*codim=*/0);
2343
2344 pvt_region_idx_.resize(numCells);
2345 for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
2346 pvt_region_idx_[cellIdx] =
2347 eclProblem.pvtRegionIndex(cellIdx);
2348 }
2349 }
2350
2351 // The number of components in the model.
2352 template<typename TypeTag>
2353 int
2354 BlackoilWellModel<TypeTag>::numComponents() const
2355 {
2356 // The numComponents here does not reflect the actual number of the components in the system.
2357 // It more or less reflects the number of mass conservation equations for the well equations.
2358 // For example, in the current formulation, we do not have the polymer conservation equation
2359 // in the well equations. As a result, for an oil-water-polymer system, this function will return 2.
2360 // In some way, it makes this function appear to be confusing from its name, and we need
2361 // to revisit/revise this function again when extending the variants of system that flow can simulate.
2362 int numComp = numPhases() < 3? numPhases(): FluidSystem::numComponents;
2363 if constexpr (has_solvent_) {
2364 numComp++;
2365 }
2366 return numComp;
2367 }
2368
2369 template<typename TypeTag>
2370 void
2371 BlackoilWellModel<TypeTag>::extractLegacyDepth_()
2372 {
2373 const auto& eclProblem = ebosSimulator_.problem();
2374 depth_.resize(local_num_cells_);
2375 for (unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
2376 depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
2377 }
2378 }
2379
2380 template<typename TypeTag>
2381 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
2382 BlackoilWellModel<TypeTag>::
2383 getWell(const std::string& well_name) const
2384 {
2385 // finding the iterator of the well in wells_ecl
2386 auto well = std::find_if(well_container_.begin(),
2387 well_container_.end(),
2388 [&well_name](const WellInterfacePtr& elem)->bool {
2389 return elem->name() == well_name;
2390 });
2391
2392 assert(well != well_container_.end());
2393
2394 return *well;
2395 }
2396
2397 template<typename TypeTag>
2398 bool
2399 BlackoilWellModel<TypeTag>::
2400 hasWell(const std::string& well_name) const
2401 {
2402 return std::any_of(well_container_.begin(), well_container_.end(),
2403 [&well_name](const WellInterfacePtr& elem) -> bool
2404 {
2405 return elem->name() == well_name;
2406 });
2407 }
2408
2409
2410
2411
2412 template <typename TypeTag>
2413 int
2414 BlackoilWellModel<TypeTag>::
2415 reportStepIndex() const
2416 {
2417 return std::max(this->ebosSimulator_.episodeIndex(), 0);
2418 }
2419
2420
2421
2422
2423
2424 template<typename TypeTag>
2425 void
2426 BlackoilWellModel<TypeTag>::
2427 calcRates(const int fipnum,
2428 const int pvtreg,
2429 const std::vector<double>& production_rates,
2430 std::vector<double>& resv_coeff)
2431 {
2432 rateConverter_->calcCoeff(fipnum, pvtreg, production_rates, resv_coeff);
2433 }
2434
2435 template<typename TypeTag>
2436 void
2437 BlackoilWellModel<TypeTag>::
2438 calcInjRates(const int fipnum,
2439 const int pvtreg,
2440 std::vector<double>& resv_coeff)
2441 {
2442 rateConverter_->calcInjCoeff(fipnum, pvtreg, resv_coeff);
2443 }
2444
2445
2446 template <typename TypeTag>
2447 void
2448 BlackoilWellModel<TypeTag>::
2449 computeWellTemperature()
2450 {
2451 if (!has_energy_)
2452 return;
2453
2454 int np = numPhases();
2455 double cellInternalEnergy;
2456 double cellBinv;
2457 double cellDensity;
2458 double perfPhaseRate;
2459 const int nw = numLocalWells();
2460 for (auto wellID = 0*nw; wellID < nw; ++wellID) {
2461 const Well& well = wells_ecl_[wellID];
2462 if (well.isInjector())
2463 continue;
2464
2465 std::array<double,2> weighted{0.0,0.0};
2466 auto& [weighted_temperature, total_weight] = weighted;
2467
2468 auto& well_info = local_parallel_well_info_[wellID].get();
2469 auto& ws = this->wellState().well(wellID);
2470 auto& perf_data = ws.perf_data;
2471 auto& perf_phase_rate = perf_data.phase_rates;
2472
2473 using int_type = decltype(well_perf_data_[wellID].size());
2474 for (int_type perf = 0, end_perf = well_perf_data_[wellID].size(); perf < end_perf; ++perf) {
2475 const int cell_idx = well_perf_data_[wellID][perf].cell_index;
2476 const auto& intQuants = ebosSimulator_.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
2477 const auto& fs = intQuants.fluidState();
2478
2479 // we on only have one temperature pr cell any phaseIdx will do
2480 double cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
2481
2482 double weight_factor = 0.0;
2483 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
2484 {
2485 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2486 continue;
2487 }
2488 cellInternalEnergy = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
2489 cellBinv = fs.invB(phaseIdx).value();
2490 cellDensity = fs.density(phaseIdx).value();
2491 perfPhaseRate = perf_phase_rate[ perf*np + phaseIdx ];
2492 weight_factor += cellDensity * perfPhaseRate/cellBinv * cellInternalEnergy/cellTemperatures;
2493 }
2494 total_weight += weight_factor;
2495 weighted_temperature += weight_factor * cellTemperatures;
2496 }
2497 well_info.communication().sum(weighted.data(), 2);
2498 this->wellState().well(wellID).temperature = weighted_temperature/total_weight;
2499 }
2500 }
2501
2502
2503 template <typename TypeTag>
2504 void
2505 BlackoilWellModel<TypeTag>::
2506 logPrimaryVars() const
2507 {
2508 std::ostringstream os;
2509 for (const auto& w : well_container_) {
2510 os << w->name() << ":";
2511 auto pv = w->getPrimaryVars();
2512 for (const double v : pv) {
2513 os << ' ' << v;
2514 }
2515 os << '\n';
2516 }
2517 OpmLog::debug(os.str());
2518 }
2519
2520
2521
2522 template <typename TypeTag>
2523 std::vector<double>
2524 BlackoilWellModel<TypeTag>::
2525 getPrimaryVarsDomain(const Domain& domain) const
2526 {
2527 std::vector<double> ret;
2528 for (const auto& well : well_container_) {
2529 if (well_domain_.at(well->name()) == domain.index) {
2530 const auto& pv = well->getPrimaryVars();
2531 ret.insert(ret.end(), pv.begin(), pv.end());
2532 }
2533 }
2534 return ret;
2535 }
2536
2537
2538
2539 template <typename TypeTag>
2540 void
2541 BlackoilWellModel<TypeTag>::
2542 setPrimaryVarsDomain(const Domain& domain, const std::vector<double>& vars)
2543 {
2544 std::size_t offset = 0;
2545 for (auto& well : well_container_) {
2546 if (well_domain_.at(well->name()) == domain.index) {
2547 int num_pri_vars = well->setPrimaryVars(vars.begin() + offset);
2548 offset += num_pri_vars;
2549 }
2550 }
2551 assert(offset == vars.size());
2552 }
2553
2554
2555
2556 template <typename TypeTag>
2557 void
2558 BlackoilWellModel<TypeTag>::
2559 setupDomains(const std::vector<Domain>& domains)
2560 {
2561 OPM_BEGIN_PARALLEL_TRY_CATCH();
2562 // TODO: This loop nest may be slow for very large numbers
2563 // of domains and wells, but that has not been observed on
2564 // tests so far. Using the partition vector instead would
2565 // be faster if we need to change.
2566 for (const auto& wellPtr : this->well_container_) {
2567 const int first_well_cell = wellPtr->cells()[0];
2568 for (const auto& domain : domains) {
2569 const bool found = std::binary_search(domain.cells.begin(), domain.cells.end(), first_well_cell);
2570 if (found) {
2571 // Assuming that if the first well cell is found in a domain,
2572 // then all of that well's cells are in that same domain.
2573 well_domain_[wellPtr->name()] = domain.index;
2574 // Verify that all of that well's cells are in that same domain.
2575 for (int well_cell : wellPtr->cells()) {
2576 if (!std::binary_search(domain.cells.begin(), domain.cells.end(), well_cell)) {
2577 OPM_THROW(std::runtime_error, "Well found on multiple domains.");
2578 }
2579 }
2580 }
2581 }
2582 }
2583 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::setupDomains(): well found on multiple domains.",
2584 ebosSimulator_.gridView().comm());
2585
2586 // Write well/domain info to the DBG file.
2587 if (ebosSimulator_.gridView().comm().rank() == 0) {
2588 std::ostringstream os;
2589 os << "Well name Domain\n";
2590 for (const auto& [wname, domain] : well_domain_) {
2591 os << wname << std::setw(21 - wname.size()) << domain << '\n';
2592 }
2593 OpmLog::debug(os.str());
2594 }
2595 }
2596
2597} // namespace Opm
Definition AquiferInterface.hpp:35
Class for handling the blackoil well model.
Definition BlackoilWellModel.hpp:97
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Calculating the explict quantities used in the well calculation.
Definition BlackoilWellModel_impl.hpp:1821
Definition DeferredLogger.hpp:57
void logMessages()
Log all messages to the OpmLog backends, and clear the message container.
Definition DeferredLogger.cpp:85
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:60
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
ConvergenceReport gatherConvergenceReport(const ConvergenceReport &local_report, Parallel::Communication mpi_communicator)
Create a global convergence report combining local (per-process) reports.
Definition gatherConvergenceReport.cpp:249
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition phaseUsageFromDeck.cpp:145