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