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#ifndef OPM_BLACKOILWELLMODEL_IMPL_HEADER_INCLUDED
24#define OPM_BLACKOILWELLMODEL_IMPL_HEADER_INCLUDED
25
26// Improve IDE experience
27#ifndef OPM_BLACKOILWELLMODEL_HEADER_INCLUDED
28#include <config.h>
29#include <opm/simulators/wells/BlackoilWellModel.hpp>
30#endif
31
32#include <opm/grid/utility/cartesianToCompressed.hpp>
33#include <opm/common/utility/numeric/RootFinders.hpp>
34
35#include <opm/input/eclipse/Schedule/Network/Balance.hpp>
36#include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
37#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
38#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
39#include <opm/input/eclipse/Schedule/Well/WellTestConfig.hpp>
40
41#include <opm/input/eclipse/Units/UnitSystem.hpp>
42
43#include <opm/simulators/wells/BlackoilWellModelConstraints.hpp>
44#include <opm/simulators/wells/BlackoilWellModelNldd.hpp>
45#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
46#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
47#include <opm/simulators/wells/VFPProperties.hpp>
48#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
49#include <opm/simulators/wells/WellGroupControls.hpp>
50#include <opm/simulators/wells/WellGroupHelpers.hpp>
51#include <opm/simulators/wells/TargetCalculator.hpp>
52
53#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
54#include <opm/simulators/utils/MPIPacker.hpp>
55#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
56
57#if COMPILE_GPU_BRIDGE
58#include <opm/simulators/linalg/gpubridge/WellContributions.hpp>
59#endif
60
61#include <algorithm>
62#include <cassert>
63#include <iomanip>
64#include <utility>
65#include <optional>
66
67#include <fmt/format.h>
68
69namespace Opm {
70 template<typename TypeTag>
71 BlackoilWellModel<TypeTag>::
72 BlackoilWellModel(Simulator& simulator, const PhaseUsage& phase_usage)
73 : WellConnectionModule(*this, simulator.gridView().comm())
74 , BlackoilWellModelGeneric<Scalar>(simulator.vanguard().schedule(),
75 gaslift_,
76 simulator.vanguard().summaryState(),
77 simulator.vanguard().eclState(),
79 simulator.gridView().comm())
80 , simulator_(simulator)
81 , gaslift_(this->terminal_output_, this->phase_usage_)
82 {
83 local_num_cells_ = simulator_.gridView().size(0);
84
85 // Number of cells the global grid view
86 global_num_cells_ = simulator_.vanguard().globalNumCells();
87
88 {
89 auto& parallel_wells = simulator.vanguard().parallelWells();
90
91 this->parallel_well_info_.reserve(parallel_wells.size());
92 for( const auto& name_bool : parallel_wells) {
93 this->parallel_well_info_.emplace_back(name_bool, grid().comm());
94 }
95 }
96
97 this->alternative_well_rate_init_ =
98 Parameters::Get<Parameters::AlternativeWellRateInit>();
99
100 using SourceDataSpan =
101 typename PAvgDynamicSourceData<Scalar>::template SourceDataSpan<Scalar>;
102
103 this->wbp_.initializeSources(
104 [this](const std::size_t globalIndex)
105 { return this->compressedIndexForInterior(globalIndex); },
106 [this](const int localCell, SourceDataSpan sourceTerms)
107 {
108 using Item = typename SourceDataSpan::Item;
109
110 const auto* intQuants = this->simulator_.model()
111 .cachedIntensiveQuantities(localCell, /*timeIndex = */0);
112 const auto& fs = intQuants->fluidState();
113
114 sourceTerms
115 .set(Item::PoreVol, intQuants->porosity().value() *
116 this->simulator_.model().dofTotalVolume(localCell))
117 .set(Item::Depth, this->depth_[localCell]);
118
119 constexpr auto io = FluidSystem::oilPhaseIdx;
120 constexpr auto ig = FluidSystem::gasPhaseIdx;
121 constexpr auto iw = FluidSystem::waterPhaseIdx;
122
123 // Ideally, these would be 'constexpr'.
124 const auto haveOil = FluidSystem::phaseIsActive(io);
125 const auto haveGas = FluidSystem::phaseIsActive(ig);
126 const auto haveWat = FluidSystem::phaseIsActive(iw);
127
128 auto weightedPhaseDensity = [&fs](const auto ip)
129 {
130 return fs.saturation(ip).value() * fs.density(ip).value();
131 };
132
133 if (haveOil) { sourceTerms.set(Item::Pressure, fs.pressure(io).value()); }
134 else if (haveGas) { sourceTerms.set(Item::Pressure, fs.pressure(ig).value()); }
135 else { sourceTerms.set(Item::Pressure, fs.pressure(iw).value()); }
136
137 // Strictly speaking, assumes SUM(s[p]) == 1.
138 auto rho = 0.0;
139 if (haveOil) { rho += weightedPhaseDensity(io); }
140 if (haveGas) { rho += weightedPhaseDensity(ig); }
141 if (haveWat) { rho += weightedPhaseDensity(iw); }
142
143 sourceTerms.set(Item::MixtureDensity, rho);
144 }
145 );
146 }
147
148 template<typename TypeTag>
149 BlackoilWellModel<TypeTag>::
150 BlackoilWellModel(Simulator& simulator) :
151 BlackoilWellModel(simulator, phaseUsageFromDeck(simulator.vanguard().eclState()))
152 {}
153
154
155 template<typename TypeTag>
156 void
157 BlackoilWellModel<TypeTag>::
158 init()
159 {
160 extractLegacyCellPvtRegionIndex_();
161 extractLegacyDepth_();
162
163 gravity_ = simulator_.problem().gravity()[2];
164
165 this->initial_step_ = true;
166
167 // add the eWoms auxiliary module for the wells to the list
168 simulator_.model().addAuxiliaryModule(this);
169
170 is_cell_perforated_.resize(local_num_cells_, false);
171 }
172
173
174 template<typename TypeTag>
175 void
176 BlackoilWellModel<TypeTag>::
177 initWellContainer(const int reportStepIdx)
178 {
179 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
180 + ScheduleEvents::NEW_WELL;
181 const auto& events = this->schedule()[reportStepIdx].wellgroup_events();
182 for (auto& wellPtr : this->well_container_) {
183 const bool well_opened_this_step = this->report_step_starts_ &&
184 events.hasEvent(wellPtr->name(),
185 effective_events_mask);
186 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
187 this->B_avg_, well_opened_this_step);
188 }
189 }
190
191 template<typename TypeTag>
192 void
193 BlackoilWellModel<TypeTag>::
194 beginReportStep(const int timeStepIdx)
195 {
196 DeferredLogger local_deferredLogger{};
197
198 this->report_step_starts_ = true;
199 this->report_step_start_events_ = this->schedule()[timeStepIdx].wellgroup_events();
200
201 this->rateConverter_ = std::make_unique<RateConverterType>
202 (this->phase_usage_, std::vector<int>(this->local_num_cells_, 0));
203
204 {
205 // WELPI scaling runs at start of report step.
206 const auto enableWellPIScaling = true;
207 this->initializeLocalWellStructure(timeStepIdx, enableWellPIScaling);
208 }
209
210 this->initializeGroupStructure(timeStepIdx);
211
212 const auto& comm = this->simulator_.vanguard().grid().comm();
213
214 OPM_BEGIN_PARALLEL_TRY_CATCH()
215 {
216 // Create facility for calculating reservoir voidage volumes for
217 // purpose of RESV controls.
218 this->rateConverter_->template defineState<ElementContext>(this->simulator_);
219
220 // Update VFP properties.
221 {
222 const auto& sched_state = this->schedule()[timeStepIdx];
223
224 this->vfp_properties_ = std::make_unique<VFPProperties<Scalar>>
225 (sched_state.vfpinj(), sched_state.vfpprod(), this->wellState());
226 }
227 }
228 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
229 "beginReportStep() failed: ",
230 this->terminal_output_, comm)
231
232 // Store the current well and group states in order to recover in
233 // the case of failed iterations
234 this->commitWGState();
235
236 this->wellStructureChangedDynamically_ = false;
237 }
238
239
240
241
242
243 template <typename TypeTag>
244 void
246 initializeLocalWellStructure(const int reportStepIdx,
247 const bool enableWellPIScaling)
248 {
249 DeferredLogger local_deferredLogger{};
250
251 const auto& comm = this->simulator_.vanguard().grid().comm();
252
253 // Wells_ecl_ holds this rank's wells, both open and stopped/shut.
254 this->wells_ecl_ = this->getLocalWells(reportStepIdx);
255 this->local_parallel_well_info_ =
256 this->createLocalParallelWellInfo(this->wells_ecl_);
257
258 // At least initializeWellState() might be throw an exception in
259 // UniformTabulated2DFunction. Playing it safe by extending the
260 // scope a bit.
261 OPM_BEGIN_PARALLEL_TRY_CATCH()
262 {
263 this->initializeWellPerfData();
264 this->initializeWellState(reportStepIdx);
265 this->wbp_.initializeWBPCalculationService();
266
267 if (this->param_.use_multisegment_well_ && this->anyMSWellOpenLocal()) {
268 this->wellState().initWellStateMSWell(this->wells_ecl_, &this->prevWellState());
269 }
270
271 this->initializeWellProdIndCalculators();
272
273 if (enableWellPIScaling && this->schedule()[reportStepIdx].events()
274 .hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX))
275 {
276 this->runWellPIScaling(reportStepIdx, local_deferredLogger);
277 }
278 }
279 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
280 "Failed to initialize local well structure: ",
281 this->terminal_output_, comm)
282 }
283
284
285
286
287
288 template <typename TypeTag>
289 void
291 initializeGroupStructure(const int reportStepIdx)
292 {
293 DeferredLogger local_deferredLogger{};
294
295 const auto& comm = this->simulator_.vanguard().grid().comm();
296
297 OPM_BEGIN_PARALLEL_TRY_CATCH()
298 {
299 const auto& fieldGroup =
300 this->schedule().getGroup("FIELD", reportStepIdx);
301
303 this->schedule(),
304 this->summaryState(),
305 reportStepIdx,
306 this->groupState());
307
308 // Define per region average pressure calculators for use by
309 // pressure maintenance groups (GPMAINT keyword).
310 if (this->schedule()[reportStepIdx].has_gpmaint()) {
312 (fieldGroup,
313 this->schedule(),
314 reportStepIdx,
315 this->eclState_.fieldProps(),
316 this->phase_usage_,
317 this->regionalAveragePressureCalculator_);
318 }
319 }
320 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
321 "Failed to initialize group structure: ",
322 this->terminal_output_, comm)
323 }
324
325
326
327
328
329 // called at the beginning of a time step
330 template<typename TypeTag>
331 void
334 {
335 OPM_TIMEBLOCK(beginTimeStep);
336
337 this->updateAverageFormationFactor();
338
339 DeferredLogger local_deferredLogger;
340
341 this->switched_prod_groups_.clear();
342 this->switched_inj_groups_.clear();
343
344 if (this->wellStructureChangedDynamically_) {
345 // Something altered the well structure/topology. Possibly
346 // WELSPECS/COMPDAT and/or WELOPEN run from an ACTIONX block.
347 // Reconstruct the local wells to account for the new well
348 // structure.
349 const auto reportStepIdx =
350 this->simulator_.episodeIndex();
351
352 // Disable WELPI scaling when well structure is updated in the
353 // middle of a report step.
354 const auto enableWellPIScaling = false;
355
356 this->initializeLocalWellStructure(reportStepIdx, enableWellPIScaling);
357 this->initializeGroupStructure(reportStepIdx);
358
359 this->commitWGState();
360
361 // Reset topology flag to signal that we've handled this
362 // structure change. That way we don't end up here in
363 // subsequent calls to beginTimeStep() unless there's a new
364 // dynamic change to the well structure during a report step.
365 this->wellStructureChangedDynamically_ = false;
366 }
367
368 this->resetWGState();
369
370 const int reportStepIdx = simulator_.episodeIndex();
371 this->updateAndCommunicateGroupData(reportStepIdx,
372 simulator_.model().newtonMethod().numIterations(),
373 param_.nupcol_group_rate_tolerance_,
374 local_deferredLogger);
375
376 this->wellState().updateWellsDefaultALQ(this->schedule(), reportStepIdx, this->summaryState());
377 this->wellState().gliftTimeStepInit();
378
379 const double simulationTime = simulator_.time();
380 OPM_BEGIN_PARALLEL_TRY_CATCH();
381 {
382 // test wells
383 wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
384
385 // create the well container
386 createWellContainer(reportStepIdx);
387
388 // Wells are active if they are active wells on at least one process.
389 const Grid& grid = simulator_.vanguard().grid();
390 this->wells_active_ = grid.comm().max(!this->well_container_.empty());
391
392 // do the initialization for all the wells
393 // TODO: to see whether we can postpone of the intialization of the well containers to
394 // optimize the usage of the following several member variables
395 this->initWellContainer(reportStepIdx);
396
397 // update the updated cell flag
398 std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
399 for (auto& well : well_container_) {
400 well->updatePerforatedCell(is_cell_perforated_);
401 }
402
403 // calculate the efficiency factors for each well
404 this->calculateEfficiencyFactors(reportStepIdx);
405
406 if constexpr (has_polymer_)
407 {
408 if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
409 this->setRepRadiusPerfLength();
410 }
411 }
412
413 }
414 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "beginTimeStep() failed: ",
415 this->terminal_output_, simulator_.vanguard().grid().comm());
416
417 for (auto& well : well_container_) {
418 well->setVFPProperties(this->vfp_properties_.get());
419 well->setGuideRate(&this->guideRate_);
420 }
421
422 this->updateFiltrationModelsPreStep(local_deferredLogger);
423
424 // Close completions due to economic reasons
425 for (auto& well : well_container_) {
426 well->closeCompletions(this->wellTestState());
427 }
428
429 // we need the inj_multiplier from the previous time step
430 this->initInjMult();
431
432 const auto& summaryState = simulator_.vanguard().summaryState();
433 if (alternative_well_rate_init_) {
434 // Update the well rates of well_state_, if only single-phase rates, to
435 // have proper multi-phase rates proportional to rates at bhp zero.
436 // This is done only for producers, as injectors will only have a single
437 // nonzero phase anyway.
438 for (const auto& well : well_container_) {
439 const bool zero_target = well->stoppedOrZeroRateTarget(simulator_, this->wellState(), local_deferredLogger);
440 if (well->isProducer() && !zero_target) {
441 well->updateWellStateRates(simulator_, this->wellState(), local_deferredLogger);
442 }
443 }
444 }
445
446 for (const auto& well : well_container_) {
447 if (well->isVFPActive(local_deferredLogger)){
448 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
449 }
450 }
451
452 // calculate the well potentials
453 try {
454 this->updateWellPotentials(reportStepIdx,
455 /*onlyAfterEvent*/true,
456 simulator_.vanguard().summaryConfig(),
457 local_deferredLogger);
458 } catch ( std::runtime_error& e ) {
459 const std::string msg = "A zero well potential is returned for output purposes. ";
460 local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
461 }
462
463 //update guide rates
464 const auto& comm = simulator_.vanguard().grid().comm();
465 std::vector<Scalar> pot(this->numPhases(), 0.0);
466 const Group& fieldGroup = this->schedule().getGroup("FIELD", reportStepIdx);
467 WellGroupHelpers<Scalar>::updateGuideRates(fieldGroup,
468 this->schedule(),
469 summaryState,
470 this->phase_usage_,
471 reportStepIdx,
472 simulationTime,
473 this->wellState(),
474 this->groupState(),
475 comm,
476 &this->guideRate_,
477 pot,
478 local_deferredLogger);
479 std::string exc_msg;
480 auto exc_type = ExceptionType::NONE;
481 // update gpmaint targets
482 if (this->schedule_[reportStepIdx].has_gpmaint()) {
483 for (const auto& calculator : regionalAveragePressureCalculator_) {
484 calculator.second->template defineState<ElementContext>(simulator_);
485 }
486 const double dt = simulator_.timeStepSize();
487 WellGroupHelpers<Scalar>::updateGpMaintTargetForGroups(fieldGroup,
488 this->schedule_,
489 regionalAveragePressureCalculator_,
490 reportStepIdx,
491 dt,
492 this->wellState(),
493 this->groupState());
494 }
495 try {
496 // Compute initial well solution for new wells and injectors that change injection type i.e. WAG.
497 for (auto& well : well_container_) {
498 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
499 + ScheduleEvents::INJECTION_TYPE_CHANGED
500 + ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
501 + ScheduleEvents::NEW_WELL;
502
503 const auto& events = this->schedule()[reportStepIdx].wellgroup_events();
504 const bool event = this->report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
505 const bool dyn_status_change = this->wellState().well(well->name()).status
506 != this->prevWellState().well(well->name()).status;
507
508 if (event || dyn_status_change) {
509 try {
510 well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), local_deferredLogger);
511 well->calculateExplicitQuantities(simulator_, this->wellState(), local_deferredLogger);
512 well->solveWellEquation(simulator_, this->wellState(), this->groupState(), local_deferredLogger);
513 } catch (const std::exception& e) {
514 const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates";
515 local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
516 }
517 }
518 }
519 }
520 // Catch clauses for all errors setting exc_type and exc_msg
521 OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
522
523 if (exc_type != ExceptionType::NONE) {
524 const std::string msg = "Compute initial well solution for new wells failed. Continue with zero initial rates";
525 local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
526 }
527
528 logAndCheckForExceptionsAndThrow(local_deferredLogger,
529 exc_type, "beginTimeStep() failed: " + exc_msg, this->terminal_output_, comm);
530
531 }
532
533 template<typename TypeTag>
534 void
535 BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx,
536 const double simulationTime,
537 DeferredLogger& deferred_logger)
538 {
539 for (const std::string& well_name : this->getWellsForTesting(timeStepIdx, simulationTime)) {
540 const Well& wellEcl = this->schedule().getWell(well_name, timeStepIdx);
541 if (wellEcl.getStatus() == Well::Status::SHUT)
542 continue;
543
544 WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
545 // some preparation before the well can be used
546 well->init(&this->phase_usage_, depth_, gravity_, B_avg_, true);
547
548 Scalar well_efficiency_factor = wellEcl.getEfficiencyFactor() *
549 this->wellState().getGlobalEfficiencyScalingFactor(well_name);
550 WellGroupHelpers<Scalar>::accumulateGroupEfficiencyFactor(this->schedule().getGroup(wellEcl.groupName(),
551 timeStepIdx),
552 this->schedule(),
553 timeStepIdx,
554 well_efficiency_factor);
555
556 well->setWellEfficiencyFactor(well_efficiency_factor);
557 well->setVFPProperties(this->vfp_properties_.get());
558 well->setGuideRate(&this->guideRate_);
559
560 // initialize rates/previous rates to prevent zero fractions in vfp-interpolation
561 if (well->isProducer()) {
562 well->updateWellStateRates(simulator_, this->wellState(), deferred_logger);
563 }
564 if (well->isVFPActive(deferred_logger)) {
565 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
566 }
567
568 const auto& network = this->schedule()[timeStepIdx].network();
569 if (network.active() && !this->node_pressures_.empty()) {
570 if (well->isProducer()) {
571 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
572 if (it != this->node_pressures_.end()) {
573 // The well belongs to a group which has a network nodal pressure,
574 // set the dynamic THP constraint based on the network nodal pressure
575 const Scalar nodal_pressure = it->second;
576 well->setDynamicThpLimit(nodal_pressure);
577 }
578 }
579 }
580 try {
581 using GLiftEclWells = typename GasLiftGroupInfo<Scalar>::GLiftEclWells;
582 GLiftEclWells ecl_well_map;
583 gaslift_.initGliftEclWellMap(well_container_, ecl_well_map);
584 well->wellTesting(simulator_,
585 simulationTime,
586 this->wellState(),
587 this->groupState(),
588 this->wellTestState(),
589 this->phase_usage_,
590 ecl_well_map,
591 this->well_open_times_,
592 deferred_logger);
593 } catch (const std::exception& e) {
594 const std::string msg = fmt::format("Exception during testing of well: {}. The well will not open.\n Exception message: {}", wellEcl.name(), e.what());
595 deferred_logger.warning("WELL_TESTING_FAILED", msg);
596 }
597 }
598 }
599
600
601
602
603
604 // called at the end of a report step
605 template<typename TypeTag>
606 void
607 BlackoilWellModel<TypeTag>::
608 endReportStep()
609 {
610 // Clear the communication data structures for above values.
611 for (auto&& pinfo : this->local_parallel_well_info_)
612 {
613 pinfo.get().clear();
614 }
615 }
616
617
618
619
620
621 // called at the end of a report step
622 template<typename TypeTag>
623 const SimulatorReportSingle&
624 BlackoilWellModel<TypeTag>::
625 lastReport() const {return last_report_; }
626
627
628
629
630
631 // called at the end of a time step
632 template<typename TypeTag>
633 void
634 BlackoilWellModel<TypeTag>::
635 timeStepSucceeded(const double simulationTime, const double dt)
636 {
637 this->closed_this_step_.clear();
638
639 // time step is finished and we are not any more at the beginning of an report step
640 this->report_step_starts_ = false;
641 const int reportStepIdx = simulator_.episodeIndex();
642
643 DeferredLogger local_deferredLogger;
644 for (const auto& well : well_container_) {
645 if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
646 well->updateWaterThroughput(dt, this->wellState());
647 }
648 }
649 // update connection transmissibility factor and d factor (if applicable) in the wellstate
650 for (const auto& well : well_container_) {
651 well->updateConnectionTransmissibilityFactor(simulator_, this->wellState().well(well->indexOfWell()));
652 well->updateConnectionDFactor(simulator_, this->wellState().well(well->indexOfWell()));
653 }
654
655 if (Indices::waterEnabled) {
656 this->updateFiltrationModelsPostStep(dt, FluidSystem::waterPhaseIdx, local_deferredLogger);
657 }
658
659 // WINJMULT: At the end of the time step, update the inj_multiplier saved in WellState for later use
660 this->updateInjMult(local_deferredLogger);
661
662 // report well switching
663 for (const auto& well : well_container_) {
664 well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);
665 }
666 // report group switching
667 if (this->terminal_output_) {
668 this->reportGroupSwitching(local_deferredLogger);
669 }
670
671 // update the rate converter with current averages pressures etc in
672 rateConverter_->template defineState<ElementContext>(simulator_);
673
674 // calculate the well potentials
675 try {
676 this->updateWellPotentials(reportStepIdx,
677 /*onlyAfterEvent*/false,
678 simulator_.vanguard().summaryConfig(),
679 local_deferredLogger);
680 } catch ( std::runtime_error& e ) {
681 const std::string msg = "A zero well potential is returned for output purposes. ";
682 local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
683 }
684
685 updateWellTestState(simulationTime, this->wellTestState());
686
687 // check group sales limits at the end of the timestep
688 const Group& fieldGroup = this->schedule_.getGroup("FIELD", reportStepIdx);
689 this->checkGEconLimits(fieldGroup, simulationTime,
690 simulator_.episodeIndex(), local_deferredLogger);
691 this->checkGconsaleLimits(fieldGroup, this->wellState(),
692 simulator_.episodeIndex(), local_deferredLogger);
693
694 this->calculateProductivityIndexValues(local_deferredLogger);
695
696 this->commitWGState();
697
698 const Opm::Parallel::Communication& comm = grid().comm();
699 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
700 if (this->terminal_output_) {
701 global_deferredLogger.logMessages();
702 }
703
704 //reporting output temperatures
705 this->computeWellTemperature();
706 }
707
708
709 template<typename TypeTag>
710 void
711 BlackoilWellModel<TypeTag>::
712 computeTotalRatesForDof(RateVector& rate,
713 unsigned elemIdx) const
714 {
715 rate = 0;
716
717 if (!is_cell_perforated_[elemIdx]) {
718 return;
719 }
720
721 for (const auto& well : well_container_)
722 well->addCellRates(rate, elemIdx);
723 }
724
725
726 template<typename TypeTag>
727 template <class Context>
728 void
729 BlackoilWellModel<TypeTag>::
730 computeTotalRatesForDof(RateVector& rate,
731 const Context& context,
732 unsigned spaceIdx,
733 unsigned timeIdx) const
734 {
735 rate = 0;
736 int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
737
738 if (!is_cell_perforated_[elemIdx]) {
739 return;
740 }
741
742 for (const auto& well : well_container_)
743 well->addCellRates(rate, elemIdx);
744 }
745
746
747
748 template<typename TypeTag>
749 void
750 BlackoilWellModel<TypeTag>::
751 initializeWellState(const int timeStepIdx)
752 {
753 const auto pressIx = []()
754 {
755 if (Indices::oilEnabled) { return FluidSystem::oilPhaseIdx; }
756 if (Indices::waterEnabled) { return FluidSystem::waterPhaseIdx; }
757
758 return FluidSystem::gasPhaseIdx;
759 }();
760
761 auto cellPressures = std::vector<Scalar>(this->local_num_cells_, Scalar{0});
762 auto cellTemperatures = std::vector<Scalar>(this->local_num_cells_, Scalar{0});
763
764 auto elemCtx = ElementContext { this->simulator_ };
765 const auto& gridView = this->simulator_.vanguard().gridView();
766
767 OPM_BEGIN_PARALLEL_TRY_CATCH();
768 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
769 elemCtx.updatePrimaryStencil(elem);
770 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
771
772 const auto ix = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
773 const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState();
774
775 cellPressures[ix] = fs.pressure(pressIx).value();
776 cellTemperatures[ix] = fs.temperature(0).value();
777 }
778 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ",
779 this->simulator_.vanguard().grid().comm());
780
781 this->wellState().init(cellPressures, cellTemperatures, this->schedule(), this->wells_ecl_,
782 this->local_parallel_well_info_, timeStepIdx,
783 &this->prevWellState(), this->well_perf_data_,
784 this->summaryState(), simulator_.vanguard().enableDistributedWells());
785 }
786
787
788
789
790
791 template<typename TypeTag>
792 void
793 BlackoilWellModel<TypeTag>::
794 createWellContainer(const int report_step)
795 {
796 DeferredLogger local_deferredLogger;
797
798 const int nw = this->numLocalWells();
799
800 well_container_.clear();
801
802 if (nw > 0) {
803 well_container_.reserve(nw);
804
805 const auto& wmatcher = this->schedule().wellMatcher(report_step);
806 const auto& wcycle = this->schedule()[report_step].wcycle.get();
807
808 // First loop and check for status changes. This is necessary
809 // as wcycle needs the updated open/close times.
810 std::for_each(this->wells_ecl_.begin(), this->wells_ecl_.end(),
811 [this, &wg_events = this->report_step_start_events_](const auto& well_ecl)
812 {
813 if (!well_ecl.hasConnections()) {
814 // No connections in this well. Nothing to do.
815 return;
816 }
817
818 constexpr auto events_mask = ScheduleEvents::WELL_STATUS_CHANGE |
819 ScheduleEvents::REQUEST_OPEN_WELL;
820 const bool well_status_change =
821 this->report_step_starts_ &&
822 wg_events.hasEvent(well_ecl.name(), events_mask);
823 if (well_status_change) {
824 if (well_ecl.getStatus() == WellStatus::OPEN) {
825 this->well_open_times_.insert_or_assign(well_ecl.name(),
826 this->simulator_.time());
827 this->well_close_times_.erase(well_ecl.name());
828 } else if (well_ecl.getStatus() == WellStatus::SHUT) {
829 this->well_close_times_.insert_or_assign(well_ecl.name(),
830 this->simulator_.time());
831 this->well_open_times_.erase(well_ecl.name());
832 }
833 }
834 });
835
836 // Grab wcycle states. This needs to run before the schedule gets processed
837 const auto cycle_states = wcycle.wellStatus(this->simulator_.time(),
838 wmatcher,
839 this->well_open_times_,
840 this->well_close_times_);
841
842 for (int w = 0; w < nw; ++w) {
843 const Well& well_ecl = this->wells_ecl_[w];
844
845 if (!well_ecl.hasConnections()) {
846 // No connections in this well. Nothing to do.
847 continue;
848 }
849
850 const std::string& well_name = well_ecl.name();
851 const auto well_status = this->schedule()
852 .getWell(well_name, report_step).getStatus();
853
854 if ((well_ecl.getStatus() == Well::Status::SHUT) ||
855 (well_status == Well::Status::SHUT))
856 {
857 // Due to ACTIONX the well might have been closed behind our back.
858 if (well_ecl.getStatus() != Well::Status::SHUT) {
859 this->closed_this_step_.insert(well_name);
860 this->wellState().shutWell(w);
861 }
862
863 this->well_open_times_.erase(well_name);
864 this->well_close_times_.erase(well_name);
865 continue;
866 }
867
868 // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
869 if (this->wellTestState().well_is_closed(well_name)) {
870 // The well was shut this timestep, we are most likely retrying
871 // a timestep without the well in question, after it caused
872 // repeated timestep cuts. It should therefore not be opened,
873 // even if it was new or received new targets this report step.
874 const bool closed_this_step = (this->wellTestState().lastTestTime(well_name) == simulator_.time());
875 // TODO: more checking here, to make sure this standard more specific and complete
876 // maybe there is some WCON keywords will not open the well
877 auto& events = this->wellState().well(w).events;
878 if (events.hasEvent(ScheduleEvents::REQUEST_OPEN_WELL)) {
879 if (!closed_this_step) {
880 this->wellTestState().open_well(well_name);
881 this->wellTestState().open_completions(well_name);
882 this->well_open_times_.insert_or_assign(well_name,
883 this->simulator_.time());
884 this->well_close_times_.erase(well_name);
885 }
886 events.clearEvent(ScheduleEvents::REQUEST_OPEN_WELL);
887 }
888 }
889
890 // TODO: should we do this for all kinds of closing reasons?
891 // something like wellTestState().hasWell(well_name)?
892 bool wellIsStopped = false;
893 if (this->wellTestState().well_is_closed(well_name))
894 {
895 if (well_ecl.getAutomaticShutIn()) {
896 // shut wells are not added to the well container
897 this->wellState().shutWell(w);
898 this->well_close_times_.erase(well_name);
899 this->well_open_times_.erase(well_name);
900 continue;
901 } else {
902 if (!well_ecl.getAllowCrossFlow()) {
903 // stopped wells where cross flow is not allowed
904 // are not added to the well container
905 this->wellState().shutWell(w);
906 this->well_close_times_.erase(well_name);
907 this->well_open_times_.erase(well_name);
908 continue;
909 }
910 // stopped wells are added to the container but marked as stopped
911 this->wellState().stopWell(w);
912 wellIsStopped = true;
913 }
914 }
915
916 // shut wells with zero rante constraints and disallowing
917 if (!well_ecl.getAllowCrossFlow()) {
918 const bool any_zero_rate_constraint = well_ecl.isProducer()
919 ? well_ecl.productionControls(this->summaryState_).anyZeroRateConstraint()
920 : well_ecl.injectionControls(this->summaryState_).anyZeroRateConstraint();
921 if (any_zero_rate_constraint) {
922 // Treat as shut, do not add to container.
923 local_deferredLogger.debug(fmt::format(" Well {} gets shut due to having zero rate constraint and disallowing crossflow ", well_ecl.name()) );
924 this->wellState().shutWell(w);
925 this->well_close_times_.erase(well_name);
926 this->well_open_times_.erase(well_name);
927 continue;
928 }
929 }
930
931 if (well_status == Well::Status::STOP) {
932 this->wellState().stopWell(w);
933 this->well_close_times_.erase(well_name);
934 this->well_open_times_.erase(well_name);
935 wellIsStopped = true;
936 }
937
938 if (!wcycle.empty()) {
939 const auto it = cycle_states.find(well_name);
940 if (it != cycle_states.end()) {
941 if (!it->second) {
942 this->wellState().shutWell(w);
943 continue;
944 } else {
945 this->wellState().openWell(w);
946 }
947 }
948 }
949
950 well_container_.emplace_back(this->createWellPointer(w, report_step));
951
952 if (wellIsStopped) {
953 well_container_.back()->stopWell();
954 this->well_close_times_.erase(well_name);
955 this->well_open_times_.erase(well_name);
956 }
957 }
958
959 if (!wcycle.empty()) {
960 const auto schedule_open =
961 [&wg_events = this->report_step_start_events_](const std::string& name)
962 {
963 return wg_events.hasEvent(name, ScheduleEvents::REQUEST_OPEN_WELL);
964 };
965 for (const auto& [wname, wscale] : wcycle.efficiencyScale(this->simulator_.time(),
966 this->simulator_.timeStepSize(),
967 wmatcher,
968 this->well_open_times_,
969 schedule_open))
970 {
971 this->wellState().updateEfficiencyScalingFactor(wname, wscale);
972 this->schedule_.add_event(ScheduleEvents::WELLGROUP_EFFICIENCY_UPDATE, report_step);
973 }
974 }
975 }
976
977 // Collect log messages and print.
978
979 const Opm::Parallel::Communication& comm = grid().comm();
980 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
981 if (this->terminal_output_) {
982 global_deferredLogger.logMessages();
983 }
984
985 this->well_container_generic_.clear();
986 for (auto& w : well_container_)
987 this->well_container_generic_.push_back(w.get());
988
989 const auto& network = this->schedule()[report_step].network();
990 if (network.active() && !this->node_pressures_.empty()) {
991 for (auto& well: this->well_container_generic_) {
992 // Producers only, since we so far only support the
993 // "extended" network model (properties defined by
994 // BRANPROP and NODEPROP) which only applies to producers.
995 if (well->isProducer()) {
996 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
997 if (it != this->node_pressures_.end()) {
998 // The well belongs to a group which has a network nodal pressure,
999 // set the dynamic THP constraint based on the network nodal pressure
1000 const Scalar nodal_pressure = it->second;
1001 well->setDynamicThpLimit(nodal_pressure);
1002 }
1003 }
1004 }
1005 }
1006
1007 this->wbp_.registerOpenWellsForWBPCalculation();
1008 }
1009
1010
1011
1012
1013
1014 template <typename TypeTag>
1015 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1016 BlackoilWellModel<TypeTag>::
1017 createWellPointer(const int wellID, const int report_step) const
1018 {
1019 const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
1020
1021 if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
1022 return this->template createTypedWellPointer<StandardWell<TypeTag>>(wellID, report_step);
1023 }
1024 else {
1025 return this->template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, report_step);
1026 }
1027 }
1028
1029
1030
1031
1032
1033 template <typename TypeTag>
1034 template <typename WellType>
1035 std::unique_ptr<WellType>
1036 BlackoilWellModel<TypeTag>::
1037 createTypedWellPointer(const int wellID, const int time_step) const
1038 {
1039 // Use the pvtRegionIdx from the top cell
1040 const auto& perf_data = this->well_perf_data_[wellID];
1041
1042 // Cater for case where local part might have no perforations.
1043 const auto pvtreg = perf_data.empty()
1044 ? 0 : this->pvt_region_idx_[perf_data.front().cell_index];
1045
1046 const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get();
1047 const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
1048
1049 return std::make_unique<WellType>(this->wells_ecl_[wellID],
1050 parallel_well_info,
1051 time_step,
1052 this->param_,
1053 *this->rateConverter_,
1054 global_pvtreg,
1055 this->numComponents(),
1056 this->numPhases(),
1057 wellID,
1058 perf_data);
1059 }
1060
1061
1062
1063
1064
1065 template<typename TypeTag>
1066 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1067 BlackoilWellModel<TypeTag>::
1068 createWellForWellTest(const std::string& well_name,
1069 const int report_step,
1070 DeferredLogger& deferred_logger) const
1071 {
1072 // Finding the location of the well in wells_ecl
1073 const auto it = std::find_if(this->wells_ecl_.begin(),
1074 this->wells_ecl_.end(),
1075 [&well_name](const auto& w)
1076 { return well_name == w.name(); });
1077 // It should be able to find in wells_ecl.
1078 if (it == this->wells_ecl_.end()) {
1079 OPM_DEFLOG_THROW(std::logic_error,
1080 fmt::format("Could not find well {} in wells_ecl ", well_name),
1081 deferred_logger);
1082 }
1083
1084 const int pos = static_cast<int>(std::distance(this->wells_ecl_.begin(), it));
1085 return this->createWellPointer(pos, report_step);
1086 }
1087
1088
1089
1090 template<typename TypeTag>
1091 void
1092 BlackoilWellModel<TypeTag>::
1093 doPreStepNetworkRebalance(DeferredLogger& deferred_logger)
1094 {
1095 OPM_TIMEFUNCTION();
1096 const double dt = this->simulator_.timeStepSize();
1097 // TODO: should we also have the group and network backed-up here in case the solution did not get converged?
1098 auto& well_state = this->wellState();
1099 const std::size_t max_iter = param_.network_max_outer_iterations_;
1100 bool converged = false;
1101 std::size_t iter = 0;
1102 bool changed_well_group = false;
1103 do {
1104 changed_well_group = updateWellControlsAndNetwork(true, dt, deferred_logger);
1105 assembleWellEqWithoutIteration(dt, deferred_logger);
1106 converged = this->getWellConvergence(this->B_avg_, true).converged() && !changed_well_group;
1107 if (converged) {
1108 break;
1109 }
1110 ++iter;
1111 OPM_BEGIN_PARALLEL_TRY_CATCH();
1112 for (auto& well : this->well_container_) {
1113 well->solveEqAndUpdateWellState(simulator_, well_state, deferred_logger);
1114 }
1115 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::doPreStepNetworkRebalance() failed: ",
1116 this->simulator_.vanguard().grid().comm());
1117 } while (iter < max_iter);
1118
1119 if (!converged) {
1120 const std::string msg = fmt::format("Initial (pre-step) network balance did not get converged with {} iterations, "
1121 "unconverged network balance result will be used", max_iter);
1122 deferred_logger.warning(msg);
1123 } else {
1124 const std::string msg = fmt::format("Initial (pre-step) network balance converged with {} iterations", iter);
1125 deferred_logger.debug(msg);
1126 }
1127 }
1128
1129
1130
1131
1132 template<typename TypeTag>
1133 void
1134 BlackoilWellModel<TypeTag>::
1135 assemble(const int iterationIdx,
1136 const double dt)
1137 {
1138 OPM_TIMEFUNCTION();
1139 DeferredLogger local_deferredLogger;
1140
1141 if constexpr (BlackoilWellModelGasLift<TypeTag>::glift_debug) {
1142 if (gaslift_.terminalOutput()) {
1143 const std::string msg =
1144 fmt::format("assemble() : iteration {}" , iterationIdx);
1145 gaslift_.gliftDebug(msg, local_deferredLogger);
1146 }
1147 }
1148 last_report_ = SimulatorReportSingle();
1149 Dune::Timer perfTimer;
1150 perfTimer.start();
1151 this->closed_offending_wells_.clear();
1152
1153 {
1154 const int episodeIdx = simulator_.episodeIndex();
1155 const auto& network = this->schedule()[episodeIdx].network();
1156 if (!this->wellsActive() && !network.active()) {
1157 return;
1158 }
1159 }
1160
1161 if (iterationIdx == 0 && this->wellsActive()) {
1162 OPM_TIMEBLOCK(firstIterationAssmble);
1163 // try-catch is needed here as updateWellControls
1164 // contains global communication and has either to
1165 // be reached by all processes or all need to abort
1166 // before.
1167 OPM_BEGIN_PARALLEL_TRY_CATCH();
1168 {
1169 calculateExplicitQuantities(local_deferredLogger);
1170 prepareTimeStep(local_deferredLogger);
1171 }
1172 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1173 "assemble() failed (It=0): ",
1174 this->terminal_output_, grid().comm());
1175 }
1176
1177 const bool well_group_control_changed = updateWellControlsAndNetwork(false, dt, local_deferredLogger);
1178
1179 // even when there is no wells active, the network nodal pressure still need to be updated through updateWellControlsAndNetwork()
1180 // but there is no need to assemble the well equations
1181 if ( ! this->wellsActive() ) {
1182 return;
1183 }
1184
1185 assembleWellEqWithoutIteration(dt, local_deferredLogger);
1186
1187 // if group or well control changes we don't consider the
1188 // case converged
1189 last_report_.well_group_control_changed = well_group_control_changed;
1190 last_report_.assemble_time_well += perfTimer.stop();
1191 }
1192
1193
1194
1195
1196 template<typename TypeTag>
1197 bool
1198 BlackoilWellModel<TypeTag>::
1199 updateWellControlsAndNetwork(const bool mandatory_network_balance,
1200 const double dt,
1201 DeferredLogger& local_deferredLogger)
1202 {
1203 OPM_TIMEFUNCTION();
1204 // not necessarily that we always need to update once of the network solutions
1205 bool do_network_update = true;
1206 bool well_group_control_changed = false;
1207 Scalar network_imbalance = 0.0;
1208 // after certain number of the iterations, we use relaxed tolerance for the network update
1209 const std::size_t iteration_to_relax = param_.network_max_strict_outer_iterations_;
1210 // after certain number of the iterations, we terminate
1211 const std::size_t max_iteration = param_.network_max_outer_iterations_;
1212 std::size_t network_update_iteration = 0;
1213 while (do_network_update) {
1214 if (network_update_iteration >= max_iteration ) {
1215 // only output to terminal if we at the last newton iterations where we try to balance the network.
1216 const int episodeIdx = simulator_.episodeIndex();
1217 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1218 if (this->shouldBalanceNetwork(episodeIdx, iterationIdx + 1)) {
1219 const std::string msg = fmt::format("Maximum of {:d} network iterations has been used and we stop the update, \n"
1220 "and try again after the next Newton iteration (imbalance = {:.2e} bar, ctrl_change = {})",
1221 max_iteration, network_imbalance*1.0e-5, well_group_control_changed);
1222 local_deferredLogger.debug(msg);
1223 } else {
1224 if (this->terminal_output_) {
1225 const std::string msg = fmt::format("Maximum of {:d} network iterations has been used and we stop the update. \n"
1226 "The simulator will continue with unconverged network results (imbalance = {:.2e} bar, ctrl_change = {})",
1227 max_iteration, network_imbalance*1.0e-5, well_group_control_changed);
1228 local_deferredLogger.info(msg);
1229 }
1230 }
1231 break;
1232 }
1233 if (this->terminal_output_ && (network_update_iteration == iteration_to_relax) ) {
1234 local_deferredLogger.debug("We begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations ");
1235 }
1236 const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
1237 // Never optimize gas lift in last iteration, to allow network convergence (unless max_iter < 2)
1238 const bool optimize_gas_lift = ( (network_update_iteration + 1) < std::max(max_iteration, static_cast<std::size_t>(2)) );
1239 std::tie(well_group_control_changed, do_network_update, network_imbalance) =
1240 updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, optimize_gas_lift, dt,local_deferredLogger);
1241 ++network_update_iteration;
1242 }
1243 return well_group_control_changed;
1244 }
1245
1246
1247
1248
1249 template<typename TypeTag>
1250 std::tuple<bool, bool, typename BlackoilWellModel<TypeTag>::Scalar>
1251 BlackoilWellModel<TypeTag>::
1252 updateWellControlsAndNetworkIteration(const bool mandatory_network_balance,
1253 const bool relax_network_tolerance,
1254 const bool optimize_gas_lift,
1255 const double dt,
1256 DeferredLogger& local_deferredLogger)
1257 {
1258 OPM_TIMEFUNCTION();
1259 auto [well_group_control_changed, more_inner_network_update, network_imbalance] =
1260 updateWellControls(mandatory_network_balance,
1261 local_deferredLogger,
1262 relax_network_tolerance);
1263
1264 bool alq_updated = false;
1265 OPM_BEGIN_PARALLEL_TRY_CATCH();
1266 {
1267 if (optimize_gas_lift) alq_updated = gaslift_.maybeDoGasLiftOptimize(simulator_,
1268 well_container_,
1269 this->wellState(),
1270 this->groupState(),
1271 local_deferredLogger);
1272
1273 prepareWellsBeforeAssembling(dt, local_deferredLogger);
1274 }
1275 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1276 "updateWellControlsAndNetworkIteration() failed: ",
1277 this->terminal_output_, grid().comm());
1278
1279 // update guide rates
1280 const int reportStepIdx = simulator_.episodeIndex();
1281 if (alq_updated || BlackoilWellModelGuideRates(*this).
1282 guideRateUpdateIsNeeded(reportStepIdx)) {
1283 const double simulationTime = simulator_.time();
1284 const auto& comm = simulator_.vanguard().grid().comm();
1285 const auto& summaryState = simulator_.vanguard().summaryState();
1286 std::vector<Scalar> pot(this->numPhases(), 0.0);
1287 const Group& fieldGroup = this->schedule().getGroup("FIELD", reportStepIdx);
1288 WellGroupHelpers<Scalar>::updateGuideRates(fieldGroup,
1289 this->schedule(),
1290 summaryState,
1291 this->phase_usage_,
1292 reportStepIdx,
1293 simulationTime,
1294 this->wellState(),
1295 this->groupState(),
1296 comm,
1297 &this->guideRate_,
1298 pot,
1299 local_deferredLogger);
1300 }
1301 // we need to re-iterate the network when the well group controls changed or gaslift/alq is changed or
1302 // the inner iterations are did not converge
1303 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1304 const bool more_network_update = this->shouldBalanceNetwork(reportStepIdx, iterationIdx) &&
1305 (more_inner_network_update || well_group_control_changed || alq_updated);
1306 return {well_group_control_changed, more_network_update, network_imbalance};
1307 }
1308
1309 // This function is to be used for well groups in an extended network that act as a subsea manifold
1310 // The wells of such group should have a common THP and total phase rate(s) obeying (if possible)
1311 // the well group constraint set by GCONPROD
1312 template <typename TypeTag>
1313 bool
1314 BlackoilWellModel<TypeTag>::
1315 computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger)
1316 {
1317 OPM_TIMEFUNCTION();
1318 const int reportStepIdx = this->simulator_.episodeIndex();
1319 const auto& network = this->schedule()[reportStepIdx].network();
1320 const auto& balance = this->schedule()[reportStepIdx].network_balance();
1321 const Scalar thp_tolerance = balance.thp_tolerance();
1322
1323 if (!network.active()) {
1324 return false;
1325 }
1326
1327 auto& well_state = this->wellState();
1328 auto& group_state = this->groupState();
1329
1330 bool well_group_thp_updated = false;
1331 for (const std::string& nodeName : network.node_names()) {
1332 const bool has_choke = network.node(nodeName).as_choke();
1333 if (has_choke) {
1334 const auto& summary_state = this->simulator_.vanguard().summaryState();
1335 const Group& group = this->schedule().getGroup(nodeName, reportStepIdx);
1336
1337 const auto pu = this->phase_usage_;
1338 //TODO: Auto choke combined with RESV control is not supported
1339 std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
1340 Scalar gratTargetFromSales = 0.0;
1341 if (group_state.has_grat_sales_target(group.name()))
1342 gratTargetFromSales = group_state.grat_sales_target(group.name());
1343
1344 const auto ctrl = group.productionControls(summary_state);
1345 auto cmode_tmp = ctrl.cmode;
1346 Scalar target_tmp{0.0};
1347 bool fld_none = false;
1348 if (cmode_tmp == Group::ProductionCMode::FLD || cmode_tmp == Group::ProductionCMode::NONE) {
1349 fld_none = true;
1350 // Target is set for an ancestor group. Target for autochoke group to be
1351 // derived via group guide rates
1352 const Scalar efficiencyFactor = 1.0;
1353 const Group& parentGroup = this->schedule().getGroup(group.parent(), reportStepIdx);
1354 auto target = WellGroupControls<Scalar>::getAutoChokeGroupProductionTargetRate(
1355 group.name(),
1356 parentGroup,
1357 well_state,
1358 group_state,
1359 this->schedule(),
1360 summary_state,
1361 resv_coeff,
1362 efficiencyFactor,
1363 reportStepIdx,
1364 pu,
1365 &this->guideRate_,
1366 local_deferredLogger);
1367 target_tmp = target.first;
1368 cmode_tmp = target.second;
1369 }
1370 const auto cmode = cmode_tmp;
1371 WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
1372 gratTargetFromSales, nodeName, group_state,
1373 group.has_gpmaint_control(cmode));
1374 if (!fld_none)
1375 {
1376 // Target is set for the autochoke group itself
1377 target_tmp = tcalc.groupTarget(ctrl, local_deferredLogger);
1378 }
1379
1380 const Scalar orig_target = target_tmp;
1381
1382 auto mismatch = [&] (auto group_thp) {
1383 Scalar group_rate(0.0);
1384 Scalar rate(0.0);
1385 for (auto& well : this->well_container_) {
1386 std::string well_name = well->name();
1387 auto& ws = well_state.well(well_name);
1388 if (group.hasWell(well_name)) {
1389 well->setDynamicThpLimit(group_thp);
1390 const Well& well_ecl = this->wells_ecl_[well->indexOfWell()];
1391 const auto inj_controls = Well::InjectionControls(0);
1392 const auto prod_controls = well_ecl.productionControls(summary_state);
1393 well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false);
1394 rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
1395 group_rate += rate;
1396 }
1397 }
1398 return (group_rate - orig_target)/orig_target;
1399 };
1400
1401 const auto upbranch = network.uptree_branch(nodeName);
1402 const auto it = this->node_pressures_.find((*upbranch).uptree_node());
1403 const Scalar nodal_pressure = it->second;
1404 Scalar well_group_thp = nodal_pressure;
1405
1406 std::optional<Scalar> autochoke_thp;
1407 if (auto iter = this->well_group_thp_calc_.find(nodeName); iter != this->well_group_thp_calc_.end()) {
1408 autochoke_thp = this->well_group_thp_calc_.at(nodeName);
1409 }
1410
1411 //Find an initial bracket
1412 std::array<Scalar, 2> range_initial;
1413 if (!autochoke_thp.has_value()){
1414 Scalar min_thp, max_thp;
1415 // Retrieve the terminal pressure of the associated root of the manifold group
1416 std::string node_name = nodeName;
1417 while (!network.node(node_name).terminal_pressure().has_value()) {
1418 auto branch = network.uptree_branch(node_name).value();
1419 node_name = branch.uptree_node();
1420 }
1421 min_thp = network.node(node_name).terminal_pressure().value();
1422 WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp);
1423 // Narrow down the bracket
1424 Scalar low1, high1;
1425 std::array<Scalar, 2> range = {Scalar{0.9}*min_thp, Scalar{1.1}*max_thp};
1426 std::optional<Scalar> appr_sol;
1427 WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger);
1428 min_thp = low1;
1429 max_thp = high1;
1430 range_initial = {min_thp, max_thp};
1431 }
1432
1433 if (!autochoke_thp.has_value() || autochoke_thp.value() > nodal_pressure) {
1434 // The bracket is based on the initial bracket or on a range based on a previous calculated group thp
1435 std::array<Scalar, 2> range = autochoke_thp.has_value() ?
1436 std::array<Scalar, 2>{Scalar{0.9} * autochoke_thp.value(),
1437 Scalar{1.1} * autochoke_thp.value()} : range_initial;
1438 Scalar low, high;
1439 std::optional<Scalar> approximate_solution;
1440 const Scalar tolerance1 = thp_tolerance;
1441 local_deferredLogger.debug("Using brute force search to bracket the group THP");
1442 const bool finding_bracket = WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger);
1443
1444 if (approximate_solution.has_value()) {
1445 autochoke_thp = *approximate_solution;
1446 local_deferredLogger.debug("Approximate group THP value found: " + std::to_string(autochoke_thp.value()));
1447 } else if (finding_bracket) {
1448 const Scalar tolerance2 = thp_tolerance;
1449 const int max_iteration_solve = 100;
1450 int iteration = 0;
1451 autochoke_thp = RegulaFalsiBisection<ThrowOnError>::
1452 solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration);
1453 local_deferredLogger.debug(" bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " +
1454 "iteration = " + std::to_string(iteration));
1455 local_deferredLogger.debug("Group THP value = " + std::to_string(autochoke_thp.value()));
1456 } else {
1457 autochoke_thp.reset();
1458 local_deferredLogger.debug("Group THP solve failed due to bracketing failure");
1459 }
1460 }
1461 if (autochoke_thp.has_value()) {
1462 well_group_thp_calc_[nodeName] = autochoke_thp.value();
1463 // Note: The node pressure of the auto-choke node is set to well_group_thp in computeNetworkPressures()
1464 // and must be larger or equal to the pressure of the uptree node of its branch.
1465 well_group_thp = std::max(autochoke_thp.value(), nodal_pressure);
1466 }
1467
1468 for (auto& well : this->well_container_) {
1469 std::string well_name = well->name();
1470
1471 if (well->isInjector() || !well->wellEcl().predictionMode())
1472 continue;
1473
1474 if (group.hasWell(well_name)) {
1475 well->setDynamicThpLimit(well_group_thp);
1476 }
1477 const auto& ws = this->wellState().well(well->indexOfWell());
1478 const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP;
1479 if (thp_is_limit) {
1480 well->prepareWellBeforeAssembling(this->simulator_, dt, this->wellState(), this->groupState(), local_deferredLogger);
1481 }
1482 }
1483
1484 // Use the group THP in computeNetworkPressures().
1485 const auto& current_well_group_thp = group_state.is_autochoke_group(nodeName) ? group_state.well_group_thp(nodeName) : 1e30;
1486 if (std::abs(current_well_group_thp - well_group_thp) > balance.pressure_tolerance()) {
1487 well_group_thp_updated = true;
1488 group_state.update_well_group_thp(nodeName, well_group_thp);
1489 }
1490 }
1491 }
1492 return well_group_thp_updated;
1493 }
1494
1495 template<typename TypeTag>
1496 void
1497 BlackoilWellModel<TypeTag>::
1498 assembleWellEq(const double dt, DeferredLogger& deferred_logger)
1499 {
1500 OPM_TIMEFUNCTION();
1501 for (auto& well : well_container_) {
1502 well->assembleWellEq(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1503 }
1504 }
1505
1506
1507 template<typename TypeTag>
1508 void
1509 BlackoilWellModel<TypeTag>::
1510 prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger)
1511 {
1512 OPM_TIMEFUNCTION();
1513 for (auto& well : well_container_) {
1514 well->prepareWellBeforeAssembling(simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1515 }
1516 }
1517
1518
1519 template<typename TypeTag>
1520 void
1521 BlackoilWellModel<TypeTag>::
1522 assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger)
1523 {
1524 OPM_TIMEFUNCTION();
1525 // We make sure that all processes throw in case there is an exception
1526 // on one of them (WetGasPvt::saturationPressure might throw if not converged)
1527 OPM_BEGIN_PARALLEL_TRY_CATCH();
1528
1529 for (auto& well: well_container_) {
1530 well->assembleWellEqWithoutIteration(simulator_, dt, this->wellState(), this->groupState(),
1531 deferred_logger);
1532 }
1533 OPM_END_PARALLEL_TRY_CATCH_LOG(deferred_logger, "BlackoilWellModel::assembleWellEqWithoutIteration failed: ",
1534 this->terminal_output_, grid().comm());
1535
1536 }
1537
1538#if COMPILE_GPU_BRIDGE
1539 template<typename TypeTag>
1540 void
1541 BlackoilWellModel<TypeTag>::
1542 getWellContributions(WellContributions<Scalar>& wellContribs) const
1543 {
1544 // prepare for StandardWells
1545 wellContribs.setBlockSize(StandardWell<TypeTag>::Indices::numEq, StandardWell<TypeTag>::numStaticWellEq);
1546
1547 for(unsigned int i = 0; i < well_container_.size(); i++){
1548 auto& well = well_container_[i];
1549 std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1550 if (derived) {
1551 wellContribs.addNumBlocks(derived->linSys().getNumBlocks());
1552 }
1553 }
1554
1555 // allocate memory for data from StandardWells
1556 wellContribs.alloc();
1557
1558 for(unsigned int i = 0; i < well_container_.size(); i++){
1559 auto& well = well_container_[i];
1560 // maybe WellInterface could implement addWellContribution()
1561 auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag>>(well);
1562 if (derived_std) {
1563 derived_std->linSys().extract(derived_std->numStaticWellEq, wellContribs);
1564 } else {
1565 auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
1566 if (derived_ms) {
1567 derived_ms->linSys().extract(wellContribs);
1568 } else {
1569 OpmLog::warning("Warning unknown type of well");
1570 }
1571 }
1572 }
1573 }
1574#endif
1575
1576 template<typename TypeTag>
1577 void
1578 BlackoilWellModel<TypeTag>::
1579 addWellContributions(SparseMatrixAdapter& jacobian) const
1580 {
1581 for ( const auto& well: well_container_ ) {
1582 well->addWellContributions(jacobian);
1583 }
1584 }
1585
1586 template<typename TypeTag>
1587 void
1588 BlackoilWellModel<TypeTag>::
1589 addWellPressureEquations(PressureMatrix& jacobian,
1590 const BVector& weights,
1591 const bool use_well_weights) const
1592 {
1593 int nw = this->numLocalWellsEnd();
1594 int rdofs = local_num_cells_;
1595 for ( int i = 0; i < nw; i++ ) {
1596 int wdof = rdofs + i;
1597 jacobian[wdof][wdof] = 1.0;// better scaling ?
1598 }
1599
1600 for (const auto& well : well_container_) {
1601 well->addWellPressureEquations(jacobian,
1602 weights,
1603 pressureVarIndex,
1604 use_well_weights,
1605 this->wellState());
1606 }
1607 }
1608
1609 template <typename TypeTag>
1610 void BlackoilWellModel<TypeTag>::
1611 addReservoirSourceTerms(GlobalEqVector& residual,
1612 const std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress) const
1613 {
1614 // NB this loop may write multiple times to the same element
1615 // if a cell is perforated by more than one well, so it should
1616 // not be OpenMP-parallelized.
1617 for (const auto& well : well_container_) {
1618 if (!well->isOperableAndSolvable() && !well->wellIsStopped()) {
1619 continue;
1620 }
1621 const auto& cells = well->cells();
1622 const auto& rates = well->connectionRates();
1623 for (unsigned perfIdx = 0; perfIdx < rates.size(); ++perfIdx) {
1624 unsigned cellIdx = cells[perfIdx];
1625 auto rate = rates[perfIdx];
1626 rate *= -1.0;
1627 VectorBlockType res(0.0);
1628 using MatrixBlockType = typename SparseMatrixAdapter::MatrixBlock;
1629 MatrixBlockType bMat(0.0);
1630 simulator_.model().linearizer().setResAndJacobi(res, bMat, rate);
1631 residual[cellIdx] += res;
1632 *diagMatAddress[cellIdx] += bMat;
1633 }
1634 }
1635 }
1636
1637
1638 template<typename TypeTag>
1639 void
1640 BlackoilWellModel<TypeTag>::
1641 addWellPressureEquationsStruct(PressureMatrix& jacobian) const
1642 {
1643 int nw = this->numLocalWellsEnd();
1644 int rdofs = local_num_cells_;
1645 for (int i = 0; i < nw; ++i) {
1646 int wdof = rdofs + i;
1647 jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
1648 }
1649 const auto wellconnections = this->getMaxWellConnections();
1650 for (int i = 0; i < nw; ++i) {
1651 const auto& perfcells = wellconnections[i];
1652 for (int perfcell : perfcells) {
1653 int wdof = rdofs + i;
1654 jacobian.entry(wdof, perfcell) = 0.0;
1655 jacobian.entry(perfcell, wdof) = 0.0;
1656 }
1657 }
1658 }
1659
1660
1661 template<typename TypeTag>
1662 void
1663 BlackoilWellModel<TypeTag>::
1664 recoverWellSolutionAndUpdateWellState(const BVector& x)
1665 {
1666 DeferredLogger local_deferredLogger;
1667 OPM_BEGIN_PARALLEL_TRY_CATCH();
1668 {
1669 for (const auto& well : well_container_) {
1670 const auto& cells = well->cells();
1671 x_local_.resize(cells.size());
1672
1673 for (size_t i = 0; i < cells.size(); ++i) {
1674 x_local_[i] = x[cells[i]];
1675 }
1676 well->recoverWellSolutionAndUpdateWellState(simulator_, x_local_,
1677 this->wellState(), local_deferredLogger);
1678 }
1679 }
1680 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1681 "recoverWellSolutionAndUpdateWellState() failed: ",
1682 this->terminal_output_, simulator_.vanguard().grid().comm());
1683 }
1684
1685
1686 template<typename TypeTag>
1687 void
1688 BlackoilWellModel<TypeTag>::
1689 recoverWellSolutionAndUpdateWellStateDomain(const BVector& x, const int domainIdx)
1690 {
1691 if (!nldd_) {
1692 OPM_THROW(std::logic_error, "Attempt to call NLDD method without a NLDD solver");
1693 }
1694
1695 return nldd_->recoverWellSolutionAndUpdateWellState(x, domainIdx);
1696 }
1697
1698
1699 template<typename TypeTag>
1700 ConvergenceReport
1701 BlackoilWellModel<TypeTag>::
1702 getWellConvergence(const std::vector<Scalar>& B_avg, bool checkWellGroupControls) const
1703 {
1704
1705 DeferredLogger local_deferredLogger;
1706 // Get global (from all processes) convergence report.
1707 ConvergenceReport local_report;
1708 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1709 for (const auto& well : well_container_) {
1710 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
1711 local_report += well->getWellConvergence(
1712 simulator_, this->wellState(), B_avg, local_deferredLogger,
1713 iterationIdx > param_.strict_outer_iter_wells_);
1714 } else {
1715 ConvergenceReport report;
1716 using CR = ConvergenceReport;
1717 report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
1718 local_report += report;
1719 }
1720 }
1721
1722 const Opm::Parallel::Communication comm = grid().comm();
1723 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
1724 ConvergenceReport report = gatherConvergenceReport(local_report, comm);
1725
1726 // the well_group_control_changed info is already communicated
1727 if (checkWellGroupControls) {
1728 report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed);
1729 }
1730
1731 if (this->terminal_output_) {
1732 global_deferredLogger.logMessages();
1733
1734 // Log debug messages for NaN or too large residuals.
1735 for (const auto& f : report.wellFailures()) {
1736 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
1737 OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1738 } else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
1739 OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1740 }
1741 }
1742 }
1743 return report;
1744 }
1745
1746
1747
1748
1749
1750 template<typename TypeTag>
1751 void
1753 calculateExplicitQuantities(DeferredLogger& deferred_logger) const
1754 {
1755 // TODO: checking isOperableAndSolvable() ?
1756 for (auto& well : well_container_) {
1757 well->calculateExplicitQuantities(simulator_, this->wellState(), deferred_logger);
1758 }
1759 }
1760
1761
1762
1763
1764
1765 template<typename TypeTag>
1766 std::tuple<bool, bool, typename BlackoilWellModel<TypeTag>::Scalar>
1768 updateWellControls(const bool mandatory_network_balance,
1769 DeferredLogger& deferred_logger,
1770 const bool relax_network_tolerance)
1771 {
1772 OPM_TIMEFUNCTION();
1773 const int episodeIdx = simulator_.episodeIndex();
1774 const auto& network = this->schedule()[episodeIdx].network();
1775 if (!this->wellsActive() && !network.active()) {
1776 return {false, false, 0.0};
1777 }
1778
1779 const int iterationIdx = simulator_.model().newtonMethod().numIterations();
1780 const auto& comm = simulator_.vanguard().grid().comm();
1781 this->updateAndCommunicateGroupData(episodeIdx, iterationIdx, param_.nupcol_group_rate_tolerance_, deferred_logger);
1782
1783 // network related
1784 Scalar network_imbalance = 0.0;
1785 bool more_network_update = false;
1786 if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
1787 OPM_TIMEBLOCK(BalanceNetowork);
1788 const double dt = this->simulator_.timeStepSize();
1789 // Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES)
1790 bool well_group_thp_updated = computeWellGroupThp(dt, deferred_logger);
1791 const int max_number_of_sub_iterations = param_.network_max_sub_iterations_;
1792 const Scalar network_pressure_update_damping_factor = param_.network_pressure_update_damping_factor_;
1793 const Scalar network_max_pressure_update = param_.network_max_pressure_update_in_bars_ * unit::barsa;
1794 bool more_network_sub_update = false;
1795 for (int i = 0; i < max_number_of_sub_iterations; i++) {
1796 const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, network_pressure_update_damping_factor, network_max_pressure_update);
1797 network_imbalance = comm.max(local_network_imbalance);
1798 const auto& balance = this->schedule()[episodeIdx].network_balance();
1799 constexpr Scalar relaxation_factor = 10.0;
1800 const Scalar tolerance = relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() : balance.pressure_tolerance();
1801 more_network_sub_update = this->networkActive() && network_imbalance > tolerance;
1802 if (!more_network_sub_update)
1803 break;
1804
1805 for (const auto& well : well_container_) {
1806 if (well->isInjector() || !well->wellEcl().predictionMode())
1807 continue;
1808
1809 const auto it = this->node_pressures_.find(well->wellEcl().groupName());
1810 if (it != this->node_pressures_.end()) {
1811 const auto& ws = this->wellState().well(well->indexOfWell());
1812 const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP;
1813 if (thp_is_limit) {
1814 well->prepareWellBeforeAssembling(this->simulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1815 }
1816 }
1817 }
1818 this->updateAndCommunicateGroupData(episodeIdx, iterationIdx, param_.nupcol_group_rate_tolerance_, deferred_logger);
1819 }
1820 more_network_update = more_network_sub_update || well_group_thp_updated;
1821 }
1822
1823 bool changed_well_group = false;
1824 // Check group individual constraints.
1825 const Group& fieldGroup = this->schedule().getGroup("FIELD", episodeIdx);
1826 changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
1827
1828 // Check wells' group constraints and communicate.
1829 bool changed_well_to_group = false;
1830 {
1831 OPM_TIMEBLOCK(UpdateWellControls);
1832 // For MS Wells a linear solve is performed below and the matrix might be singular.
1833 // We need to communicate the exception thrown to the others and rethrow.
1834 OPM_BEGIN_PARALLEL_TRY_CATCH()
1835 for (const auto& well : well_container_) {
1836 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
1837 const bool changed_well = well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1838 if (changed_well) {
1839 changed_well_to_group = changed_well || changed_well_to_group;
1840 }
1841 }
1842 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel: updating well controls failed: ",
1843 simulator_.gridView().comm());
1844 }
1845
1846 changed_well_to_group = comm.sum(static_cast<int>(changed_well_to_group));
1847 if (changed_well_to_group) {
1848 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1849 changed_well_group = true;
1850 }
1851
1852 // Check individual well constraints and communicate.
1853 bool changed_well_individual = false;
1854 {
1855 // For MS Wells a linear solve is performed below and the matrix might be singular.
1856 // We need to communicate the exception thrown to the others and rethrow.
1857 OPM_BEGIN_PARALLEL_TRY_CATCH()
1858 for (const auto& well : well_container_) {
1859 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
1860 const bool changed_well = well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1861 if (changed_well) {
1862 changed_well_individual = changed_well || changed_well_individual;
1863 }
1864 }
1865 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel: updating well controls failed: ",
1866 simulator_.gridView().comm());
1867 }
1868
1869 changed_well_individual = comm.sum(static_cast<int>(changed_well_individual));
1870 if (changed_well_individual) {
1871 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1872 changed_well_group = true;
1873 }
1874
1875 // update wsolvent fraction for REIN wells
1876 this->updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState());
1877
1878 return { changed_well_group, more_network_update, network_imbalance };
1879 }
1880
1881 template<typename TypeTag>
1882 void
1883 BlackoilWellModel<TypeTag>::
1884 updateAndCommunicate(const int reportStepIdx,
1885 const int iterationIdx,
1886 DeferredLogger& deferred_logger)
1887 {
1888 this->updateAndCommunicateGroupData(reportStepIdx,
1889 iterationIdx,
1890 param_.nupcol_group_rate_tolerance_,
1891 deferred_logger);
1892
1893 // updateWellStateWithTarget might throw for multisegment wells hence we
1894 // have a parallel try catch here to thrown on all processes.
1895 OPM_BEGIN_PARALLEL_TRY_CATCH()
1896 // if a well or group change control it affects all wells that are under the same group
1897 for (const auto& well : well_container_) {
1898 // We only want to update wells under group-control here
1899 const auto& ws = this->wellState().well(well->indexOfWell());
1900 if (ws.production_cmode == Well::ProducerCMode::GRUP ||
1901 ws.injection_cmode == Well::InjectorCMode::GRUP)
1902 {
1903 well->updateWellStateWithTarget(simulator_, this->groupState(),
1904 this->wellState(), deferred_logger);
1905 }
1906 }
1907 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updateAndCommunicate failed: ",
1908 simulator_.gridView().comm())
1909 this->updateAndCommunicateGroupData(reportStepIdx,
1910 iterationIdx,
1911 param_.nupcol_group_rate_tolerance_,
1912 deferred_logger);
1913 }
1914
1915 template<typename TypeTag>
1916 bool
1917 BlackoilWellModel<TypeTag>::
1918 updateGroupControls(const Group& group,
1919 DeferredLogger& deferred_logger,
1920 const int reportStepIdx,
1921 const int iterationIdx)
1922 {
1923 OPM_TIMEFUNCTION();
1924 bool changed = false;
1925 // restrict the number of group switches but only after nupcol iterations.
1926 const int nupcol = this->schedule()[reportStepIdx].nupcol();
1927 const int max_number_of_group_switches = iterationIdx < nupcol ? 9999 : param_.max_number_of_group_switches_;
1928 bool changed_hc = this->checkGroupHigherConstraints( group, deferred_logger, reportStepIdx, max_number_of_group_switches);
1929 if (changed_hc) {
1930 changed = true;
1931 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1932 }
1933
1934 bool changed_individual =
1935 BlackoilWellModelConstraints(*this).
1936 updateGroupIndividualControl(group,
1937 reportStepIdx,
1938 max_number_of_group_switches,
1939 this->switched_inj_groups_,
1940 this->switched_prod_groups_,
1941 this->closed_offending_wells_,
1942 this->groupState(),
1943 this->wellState(),
1944 deferred_logger);
1945
1946 if (changed_individual) {
1947 changed = true;
1948 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1949 }
1950 // call recursively down the group hierarchy
1951 for (const std::string& groupName : group.groups()) {
1952 bool changed_this = updateGroupControls( this->schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx,iterationIdx);
1953 changed = changed || changed_this;
1954 }
1955 return changed;
1956 }
1957
1958 template<typename TypeTag>
1959 void
1961 updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
1962 {
1963 OPM_TIMEFUNCTION();
1964 DeferredLogger local_deferredLogger;
1965 for (const auto& well : well_container_) {
1966 const auto& wname = well->name();
1967 const auto wasClosed = wellTestState.well_is_closed(wname);
1968 well->checkWellOperability(simulator_,
1969 this->wellState(),
1970 local_deferredLogger);
1971 const bool under_zero_target =
1972 well->wellUnderZeroGroupRateTarget(this->simulator_,
1973 this->wellState(),
1974 local_deferredLogger);
1975 well->updateWellTestState(this->wellState().well(wname),
1976 simulationTime,
1977 /*writeMessageToOPMLog=*/ true,
1978 under_zero_target,
1979 wellTestState,
1980 local_deferredLogger);
1981
1982 if (!wasClosed && wellTestState.well_is_closed(wname)) {
1983 this->closed_this_step_.insert(wname);
1984 }
1985 }
1986
1987 const Opm::Parallel::Communication comm = grid().comm();
1988 DeferredLogger global_deferredLogger =
1989 gatherDeferredLogger(local_deferredLogger, comm);
1990
1991 for (const auto& [group_name, to] : this->closed_offending_wells_) {
1992 if (this->hasOpenLocalWell(to.second) &&
1993 !this->wasDynamicallyShutThisTimeStep(to.second))
1994 {
1995 wellTestState.close_well(to.second,
1996 WellTestConfig::Reason::GROUP,
1997 simulationTime);
1998 this->updateClosedWellsThisStep(to.second);
1999 const std::string msg =
2000 fmt::format("Procedure on exceeding {} limit is WELL for group {}. "
2001 "Well {} is {}.",
2002 to.first,
2003 group_name,
2004 to.second,
2005 "shut");
2006 global_deferredLogger.info(msg);
2007 }
2008 }
2009
2010 if (this->terminal_output_) {
2011 global_deferredLogger.logMessages();
2012 }
2013 }
2014
2015
2016 template<typename TypeTag>
2017 void
2018 BlackoilWellModel<TypeTag>::computePotentials(const std::size_t widx,
2019 const WellState<Scalar>& well_state_copy,
2020 std::string& exc_msg,
2021 ExceptionType::ExcEnum& exc_type,
2022 DeferredLogger& deferred_logger)
2023 {
2024 OPM_TIMEFUNCTION();
2025 const int np = this->numPhases();
2026 std::vector<Scalar> potentials;
2027 const auto& well = well_container_[widx];
2028 std::string cur_exc_msg;
2029 auto cur_exc_type = ExceptionType::NONE;
2030 try {
2031 well->computeWellPotentials(simulator_, well_state_copy, potentials, deferred_logger);
2032 }
2033 // catch all possible exception and store type and message.
2034 OPM_PARALLEL_CATCH_CLAUSE(cur_exc_type, cur_exc_msg);
2035 if (cur_exc_type != ExceptionType::NONE) {
2036 exc_msg += fmt::format("\nFor well {}: {}", well->name(), cur_exc_msg);
2037 }
2038 exc_type = std::max(exc_type, cur_exc_type);
2039 // Store it in the well state
2040 // potentials is resized and set to zero in the beginning of well->ComputeWellPotentials
2041 // and updated only if sucessfull. i.e. the potentials are zero for exceptions
2042 auto& ws = this->wellState().well(well->indexOfWell());
2043 for (int p = 0; p < np; ++p) {
2044 // make sure the potentials are positive
2045 ws.well_potentials[p] = std::max(Scalar{0.0}, potentials[p]);
2046 }
2047 }
2048
2049
2050
2051 template <typename TypeTag>
2052 void
2053 BlackoilWellModel<TypeTag>::
2054 calculateProductivityIndexValues(DeferredLogger& deferred_logger)
2055 {
2056 for (const auto& wellPtr : this->well_container_) {
2057 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2058 }
2059 }
2060
2061
2062
2063
2064
2065 template <typename TypeTag>
2066 void
2067 BlackoilWellModel<TypeTag>::
2068 calculateProductivityIndexValuesShutWells(const int reportStepIdx,
2069 DeferredLogger& deferred_logger)
2070 {
2071 // For the purpose of computing PI/II values, it is sufficient to
2072 // construct StandardWell instances only. We don't need to form
2073 // well objects that honour the 'isMultisegment()' flag of the
2074 // corresponding "this->wells_ecl_[shutWell]".
2075
2076 for (const auto& shutWell : this->local_shut_wells_) {
2077 if (!this->wells_ecl_[shutWell].hasConnections()) {
2078 // No connections in this well. Nothing to do.
2079 continue;
2080 }
2081
2082 auto wellPtr = this->template createTypedWellPointer
2083 <StandardWell<TypeTag>>(shutWell, reportStepIdx);
2084
2085 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_, this->B_avg_, true);
2086
2087 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
2088 }
2089 }
2090
2091
2092
2093
2094
2095 template <typename TypeTag>
2096 void
2097 BlackoilWellModel<TypeTag>::
2098 calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
2099 DeferredLogger& deferred_logger)
2100 {
2101 wellPtr->updateProductivityIndex(this->simulator_,
2102 this->prod_index_calc_[wellPtr->indexOfWell()],
2103 this->wellState(),
2104 deferred_logger);
2105 }
2106
2107
2108
2109 template<typename TypeTag>
2110 void
2111 BlackoilWellModel<TypeTag>::
2112 prepareTimeStep(DeferredLogger& deferred_logger)
2113 {
2114 // Check if there is a network with active prediction wells at this time step.
2115 const auto episodeIdx = simulator_.episodeIndex();
2116 this->updateNetworkActiveState(episodeIdx);
2117
2118 // Rebalance the network initially if any wells in the network have status changes
2119 // (Need to check this before clearing events)
2120 const bool do_prestep_network_rebalance = param_.pre_solve_network_ && this->needPreStepNetworkRebalance(episodeIdx);
2121
2122 for (const auto& well : well_container_) {
2123 auto& events = this->wellState().well(well->indexOfWell()).events;
2124 if (events.hasEvent(WellState<Scalar>::event_mask)) {
2125 well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), deferred_logger);
2126 well->updatePrimaryVariables(simulator_, this->wellState(), deferred_logger);
2127 // There is no new well control change input within a report step,
2128 // so next time step, the well does not consider to have effective events anymore.
2129 events.clearEvent(WellState<Scalar>::event_mask);
2130 }
2131 // these events only work for the first time step within the report step
2132 if (events.hasEvent(ScheduleEvents::REQUEST_OPEN_WELL)) {
2133 events.clearEvent(ScheduleEvents::REQUEST_OPEN_WELL);
2134 }
2135 // solve the well equation initially to improve the initial solution of the well model
2136 if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
2137 try {
2138 well->solveWellEquation(simulator_, this->wellState(), this->groupState(), deferred_logger);
2139 } catch (const std::exception& e) {
2140 const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the previous rates";
2141 deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
2142 }
2143 }
2144 // If we're using local well solves that include control switches, they also update
2145 // operability, so reset before main iterations begin
2146 well->resetWellOperability();
2147 }
2148 updatePrimaryVariables(deferred_logger);
2149
2150 // Actually do the pre-step network rebalance, using the updated well states and initial solutions
2151 if (do_prestep_network_rebalance) doPreStepNetworkRebalance(deferred_logger);
2152 }
2153
2154 template<typename TypeTag>
2155 void
2156 BlackoilWellModel<TypeTag>::
2157 updateAverageFormationFactor()
2158 {
2159 std::vector< Scalar > B_avg(numComponents(), Scalar() );
2160 const auto& grid = simulator_.vanguard().grid();
2161 const auto& gridView = grid.leafGridView();
2162 ElementContext elemCtx(simulator_);
2163
2164 OPM_BEGIN_PARALLEL_TRY_CATCH();
2165 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
2166 elemCtx.updatePrimaryStencil(elem);
2167 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
2168
2169 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
2170 const auto& fs = intQuants.fluidState();
2171
2172 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
2173 {
2174 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2175 continue;
2176 }
2177
2178 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
2179 auto& B = B_avg[ compIdx ];
2180
2181 B += 1 / fs.invB(phaseIdx).value();
2182 }
2183 if constexpr (has_solvent_) {
2184 auto& B = B_avg[solventSaturationIdx];
2185 B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
2186 }
2187 }
2188 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updateAverageFormationFactor() failed: ", grid.comm())
2189
2190 // compute global average
2191 grid.comm().sum(B_avg.data(), B_avg.size());
2192 for (auto& bval : B_avg)
2193 {
2194 bval /= global_num_cells_;
2195 }
2196 B_avg_ = B_avg;
2197 }
2198
2199
2200
2201
2202
2203 template<typename TypeTag>
2204 void
2205 BlackoilWellModel<TypeTag>::
2206 updatePrimaryVariables(DeferredLogger& deferred_logger)
2207 {
2208 for (const auto& well : well_container_) {
2209 well->updatePrimaryVariables(simulator_, this->wellState(), deferred_logger);
2210 }
2211 }
2212
2213 template<typename TypeTag>
2214 void
2215 BlackoilWellModel<TypeTag>::extractLegacyCellPvtRegionIndex_()
2216 {
2217 const auto& grid = simulator_.vanguard().grid();
2218 const auto& eclProblem = simulator_.problem();
2219 const unsigned numCells = grid.size(/*codim=*/0);
2220
2221 this->pvt_region_idx_.resize(numCells);
2222 for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
2223 this->pvt_region_idx_[cellIdx] =
2224 eclProblem.pvtRegionIndex(cellIdx);
2225 }
2226 }
2227
2228 // The number of components in the model.
2229 template<typename TypeTag>
2230 int
2231 BlackoilWellModel<TypeTag>::numComponents() const
2232 {
2233 // The numComponents here does not reflect the actual number of the components in the system.
2234 // It more or less reflects the number of mass conservation equations for the well equations.
2235 // For example, in the current formulation, we do not have the polymer conservation equation
2236 // in the well equations. As a result, for an oil-water-polymer system, this function will return 2.
2237 // In some way, it makes this function appear to be confusing from its name, and we need
2238 // to revisit/revise this function again when extending the variants of system that flow can simulate.
2239 int numComp = this->numPhases() < 3 ? this->numPhases() : FluidSystem::numComponents;
2240 if constexpr (has_solvent_) {
2241 numComp++;
2242 }
2243 return numComp;
2244 }
2245
2246 template<typename TypeTag>
2247 void
2248 BlackoilWellModel<TypeTag>::extractLegacyDepth_()
2249 {
2250 const auto& eclProblem = simulator_.problem();
2251 depth_.resize(local_num_cells_);
2252 for (unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
2253 depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
2254 }
2255 }
2256
2257 template<typename TypeTag>
2258 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
2259 BlackoilWellModel<TypeTag>::
2260 getWell(const std::string& well_name) const
2261 {
2262 // finding the iterator of the well in wells_ecl
2263 auto well = std::find_if(well_container_.begin(),
2264 well_container_.end(),
2265 [&well_name](const WellInterfacePtr& elem)->bool {
2266 return elem->name() == well_name;
2267 });
2268
2269 assert(well != well_container_.end());
2270
2271 return *well;
2272 }
2273
2274 template <typename TypeTag>
2275 int
2276 BlackoilWellModel<TypeTag>::
2277 reportStepIndex() const
2278 {
2279 return std::max(this->simulator_.episodeIndex(), 0);
2280 }
2281
2282
2283
2284
2285
2286 template<typename TypeTag>
2287 void
2288 BlackoilWellModel<TypeTag>::
2289 calcResvCoeff(const int fipnum,
2290 const int pvtreg,
2291 const std::vector<Scalar>& production_rates,
2292 std::vector<Scalar>& resv_coeff)
2293 {
2294 rateConverter_->calcCoeff(fipnum, pvtreg, production_rates, resv_coeff);
2295 }
2296
2297 template<typename TypeTag>
2298 void
2299 BlackoilWellModel<TypeTag>::
2300 calcInjResvCoeff(const int fipnum,
2301 const int pvtreg,
2302 std::vector<Scalar>& resv_coeff)
2303 {
2304 rateConverter_->calcInjCoeff(fipnum, pvtreg, resv_coeff);
2305 }
2306
2307
2308 template <typename TypeTag>
2309 void
2310 BlackoilWellModel<TypeTag>::
2311 computeWellTemperature()
2312 {
2313 if constexpr (has_energy_) {
2314 int np = this->numPhases();
2315 Scalar cellInternalEnergy;
2316 Scalar cellBinv;
2317 Scalar cellDensity;
2318 Scalar perfPhaseRate;
2319 const int nw = this->numLocalWells();
2320 for (auto wellID = 0*nw; wellID < nw; ++wellID) {
2321 const Well& well = this->wells_ecl_[wellID];
2322 auto& ws = this->wellState().well(wellID);
2323 if (well.isInjector()) {
2324 if (ws.status != WellStatus::STOP) {
2325 this->wellState().well(wellID).temperature = well.inj_temperature();
2326 continue;
2327 }
2328 }
2329
2330 std::array<Scalar,2> weighted{0.0,0.0};
2331 auto& [weighted_temperature, total_weight] = weighted;
2332
2333 auto& well_info = this->local_parallel_well_info_[wellID].get();
2334 auto& perf_data = ws.perf_data;
2335 auto& perf_phase_rate = perf_data.phase_rates;
2336
2337 using int_type = decltype(this->well_perf_data_[wellID].size());
2338 for (int_type perf = 0, end_perf = this->well_perf_data_[wellID].size(); perf < end_perf; ++perf) {
2339 const int cell_idx = this->well_perf_data_[wellID][perf].cell_index;
2340 const auto& intQuants = simulator_.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
2341 const auto& fs = intQuants.fluidState();
2342
2343 // we on only have one temperature pr cell any phaseIdx will do
2344 Scalar cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
2345
2346 Scalar weight_factor = 0.0;
2347 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
2348 if (!FluidSystem::phaseIsActive(phaseIdx)) {
2349 continue;
2350 }
2351 cellInternalEnergy = fs.enthalpy(phaseIdx).value() -
2352 fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
2353 cellBinv = fs.invB(phaseIdx).value();
2354 cellDensity = fs.density(phaseIdx).value();
2355 perfPhaseRate = perf_phase_rate[perf*np + phaseIdx];
2356 weight_factor += cellDensity * perfPhaseRate / cellBinv * cellInternalEnergy / cellTemperatures;
2357 }
2358 weight_factor = std::abs(weight_factor) + 1e-13;
2359 total_weight += weight_factor;
2360 weighted_temperature += weight_factor * cellTemperatures;
2361 }
2362 well_info.communication().sum(weighted.data(), 2);
2363 this->wellState().well(wellID).temperature = weighted_temperature / total_weight;
2364 }
2365 }
2366 }
2367} // namespace Opm
2368
2369#endif
Class for handling the blackoil well model.
Definition BlackoilWellModel.hpp:94
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Calculating the explict quantities used in the well calculation.
Definition BlackoilWellModel_impl.hpp:1753
Definition DeferredLogger.hpp:57
void logMessages()
Log all messages to the OpmLog backends, and clear the message container.
Definition DeferredLogger.cpp:85
Definition WellGroupHelpers.hpp:50
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition WellState.hpp:66
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