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