My Project
BlackoilWellModel_impl.hpp
1/*
2 Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics.
3 Copyright 2016 - 2018 Equinor ASA.
4 Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
5 Copyright 2016 - 2018 Norce AS
6
7 This file is part of the Open Porous Media project (OPM).
8
9 OPM is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 OPM is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with OPM. If not, see <http://www.gnu.org/licenses/>.
21*/
22
23#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24#include <opm/core/props/phaseUsageFromDeck.hpp>
25#include <opm/grid/utility/cartesianToCompressed.hpp>
26#include <opm/input/eclipse/Units/UnitSystem.hpp>
27
28#include <opm/simulators/wells/VFPProperties.hpp>
29#include <opm/simulators/utils/MPIPacker.hpp>
30#include <opm/simulators/linalg/bda/WellContributions.hpp>
31
32#if HAVE_MPI
33#include <ebos/eclmpiserializer.hh>
34#endif
35
36#include <algorithm>
37#include <utility>
38
39#include <fmt/format.h>
40
41namespace Opm {
42 template<typename TypeTag>
43 BlackoilWellModel<TypeTag>::
44 BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& phase_usage)
45 : BlackoilWellModelGeneric(ebosSimulator.vanguard().schedule(),
46 ebosSimulator.vanguard().summaryState(),
47 ebosSimulator.vanguard().eclState(),
48 phase_usage,
49 ebosSimulator.gridView().comm())
50 , ebosSimulator_(ebosSimulator)
51 {
52 terminal_output_ = ((ebosSimulator.gridView().comm().rank() == 0) &&
53 EWOMS_GET_PARAM(TypeTag, bool, EnableTerminalOutput));
54
55 local_num_cells_ = ebosSimulator_.gridView().size(0);
56
57 // Number of cells the global grid view
58 global_num_cells_ = ebosSimulator_.vanguard().globalNumCells();
59
60 {
61 auto& parallel_wells = ebosSimulator.vanguard().parallelWells();
62
63 this->parallel_well_info_.reserve(parallel_wells.size());
64 for( const auto& name_bool : parallel_wells) {
65 this->parallel_well_info_.emplace_back(name_bool, grid().comm());
66 }
67 }
68
69 this->alternative_well_rate_init_ =
70 EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
71 }
72
73 template<typename TypeTag>
74 BlackoilWellModel<TypeTag>::
75 BlackoilWellModel(Simulator& ebosSimulator) :
76 BlackoilWellModel(ebosSimulator, phaseUsageFromDeck(ebosSimulator.vanguard().eclState()))
77 {}
78
79
80 template<typename TypeTag>
81 void
82 BlackoilWellModel<TypeTag>::
83 init()
84 {
85 extractLegacyCellPvtRegionIndex_();
86 extractLegacyDepth_();
87
88 gravity_ = ebosSimulator_.problem().gravity()[2];
89
90 initial_step_ = true;
91
92 // add the eWoms auxiliary module for the wells to the list
93 ebosSimulator_.model().addAuxiliaryModule(this);
94
95 is_cell_perforated_.resize(local_num_cells_, false);
96 }
97
98
99 template<typename TypeTag>
100 void
101 BlackoilWellModel<TypeTag>::
102 initWellContainer(const int reportStepIdx)
103 {
104 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
105 + ScheduleEvents::NEW_WELL;
106 const auto& events = schedule()[reportStepIdx].wellgroup_events();
107 for (auto& wellPtr : this->well_container_) {
108 const bool well_opened_this_step = report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask);
109 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
110 this->local_num_cells_, this->B_avg_, well_opened_this_step);
111 }
112 }
113
114 template<typename TypeTag>
115 void
116 BlackoilWellModel<TypeTag>::
117 addNeighbors(std::vector<NeighborSet>& neighbors) const
118 {
119 if (!param_.matrix_add_well_contributions_) {
120 return;
121 }
122
123 // Create cartesian to compressed mapping
124 const auto& schedule_wells = schedule().getWellsatEnd();
125
126 // initialize the additional cell connections introduced by wells.
127 for (const auto& well : schedule_wells)
128 {
129 std::vector<int> wellCells;
130 // All possible connections of the well
131 const auto& connectionSet = well.getConnections();
132 wellCells.reserve(connectionSet.size());
133
134 for ( size_t c=0; c < connectionSet.size(); c++ )
135 {
136 const auto& connection = connectionSet.get(c);
137 int compressed_idx = compressedIndexForInterior(connection.global_index());
138 if ( compressed_idx >= 0 ) { // Ignore connections in inactive/remote cells.
139 wellCells.push_back(compressed_idx);
140 }
141 }
142
143 for (int cellIdx : wellCells) {
144 neighbors[cellIdx].insert(wellCells.begin(),
145 wellCells.end());
146 }
147 }
148 }
149
150 template<typename TypeTag>
151 void
152 BlackoilWellModel<TypeTag>::
153 linearize(SparseMatrixAdapter& jacobian, GlobalEqVector& res)
154 {
155 if (!param_.matrix_add_well_contributions_)
156 {
157 OPM_BEGIN_PARALLEL_TRY_CATCH();
158 {
159 // if the well contributions are not supposed to be included explicitly in
160 // the matrix, we only apply the vector part of the Schur complement here.
161 for (const auto& well: well_container_) {
162 // r = r - duneC_^T * invDuneD_ * resWell_
163 well->apply(res);
164 }
165 }
166 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::linearize failed: ",
167 ebosSimulator_.gridView().comm());
168 return;
169 }
170
171 for (const auto& well: well_container_) {
172 well->addWellContributions(jacobian);
173
174 // applying the well residual to reservoir residuals
175 // r = r - duneC_^T * invDuneD_ * resWell_
176 well->apply(res);
177 }
178 }
179
180
181 template<typename TypeTag>
182 void
183 BlackoilWellModel<TypeTag>::
184 beginReportStep(const int timeStepIdx)
185 {
186 DeferredLogger local_deferredLogger;
187
188 report_step_starts_ = true;
189
190 const Grid& grid = ebosSimulator_.vanguard().grid();
191 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
192 // Make wells_ecl_ contain only this partition's wells.
193 wells_ecl_ = getLocalWells(timeStepIdx);
194 this->local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
195
196 // at least initializeWellState might be throw
197 // exception in opm-material (UniformTabulated2DFunction.hpp)
198 // playing it safe by extending the scope a bit.
199 OPM_BEGIN_PARALLEL_TRY_CATCH();
200 {
201
202 // The well state initialize bhp with the cell pressure in the top cell.
203 // We must therefore provide it with updated cell pressures
204 this->initializeWellPerfData();
205 this->initializeWellState(timeStepIdx, summaryState);
206
207 // handling MS well related
208 if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
209 this->wellState().initWellStateMSWell(wells_ecl_, &this->prevWellState());
210 }
211
212 const Group& fieldGroup = schedule().getGroup("FIELD", timeStepIdx);
213 WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState());
214
215 // Compute reservoir volumes for RESV controls.
216 rateConverter_.reset(new RateConverterType (phase_usage_,
217 std::vector<int>(local_num_cells_, 0)));
218 rateConverter_->template defineState<ElementContext>(ebosSimulator_);
219
220 // Compute regional average pressures used by gpmaint
221 if (schedule_[timeStepIdx].has_gpmaint()) {
222 const auto& fp = this->eclState_.fieldProps();
223 const auto& fipnum = fp.get_int("FIPNUM");
224 regionalAveragePressureCalculator_.reset(new AverageRegionalPressureType (phase_usage_,fipnum));
225 }
226
227 {
228 const auto& sched_state = this->schedule()[timeStepIdx];
229 // update VFP properties
230 vfp_properties_.reset(new VFPProperties( sched_state.vfpinj(), sched_state.vfpprod(), this->prevWellState()));
231 this->initializeWellProdIndCalculators();
232 if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) {
233 this->runWellPIScaling(timeStepIdx, local_deferredLogger);
234 }
235 }
236 }
237 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "beginReportStep() failed: ",
238 terminal_output_, grid.comm());
239 // Store the current well state, to be able to recover in the case of failed iterations
240 this->commitWGState();
241 }
242
243
244 // called at the beginning of a time step
245 template<typename TypeTag>
246 void
247 BlackoilWellModel<TypeTag>::
248 beginTimeStep()
249 {
250 updatePerforationIntensiveQuantities();
251 updateAverageFormationFactor();
252 DeferredLogger local_deferredLogger;
253 switched_prod_groups_.clear();
254 switched_inj_groups_.clear();
255
256 this->resetWGState();
257 const int reportStepIdx = ebosSimulator_.episodeIndex();
258 updateAndCommunicateGroupData(reportStepIdx,
259 ebosSimulator_.model().newtonMethod().numIterations());
260 this->wellState().gliftTimeStepInit();
261 const double simulationTime = ebosSimulator_.time();
262 OPM_BEGIN_PARALLEL_TRY_CATCH();
263 {
264 // test wells
265 wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
266
267 // create the well container
268 createWellContainer(reportStepIdx);
269
270 // Wells are active if they are active wells on at least one process.
271 const Grid& grid = ebosSimulator_.vanguard().grid();
272 wells_active_ = !this->well_container_.empty();
273 wells_active_ = grid.comm().max(wells_active_);
274
275 // do the initialization for all the wells
276 // TODO: to see whether we can postpone of the intialization of the well containers to
277 // optimize the usage of the following several member variables
278 this->initWellContainer(reportStepIdx);
279
280 // update the updated cell flag
281 std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
282 for (auto& well : well_container_) {
283 well->updatePerforatedCell(is_cell_perforated_);
284 }
285
286 // calculate the efficiency factors for each well
287 calculateEfficiencyFactors(reportStepIdx);
288
289 if constexpr (has_polymer_)
290 {
291 if (PolymerModule::hasPlyshlog() || getPropValue<TypeTag, Properties::EnablePolymerMW>() ) {
292 setRepRadiusPerfLength();
293 }
294 }
295
296 }
297 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "beginTimeStep() failed: ",
298 terminal_output_, ebosSimulator_.vanguard().grid().comm());
299
300 for (auto& well : well_container_) {
301 well->setVFPProperties(vfp_properties_.get());
302 well->setGuideRate(&guideRate_);
303 }
304
305 // Close completions due to economical reasons
306 for (auto& well : well_container_) {
307 well->closeCompletions(wellTestState());
308 }
309
310 if (alternative_well_rate_init_) {
311 // Update the well rates of well_state_, if only single-phase rates, to
312 // have proper multi-phase rates proportional to rates at bhp zero.
313 // This is done only for producers, as injectors will only have a single
314 // nonzero phase anyway.
315 for (auto& well : well_container_) {
316 if (well->isProducer()) {
317 well->updateWellStateRates(ebosSimulator_, this->wellState(), local_deferredLogger);
318 }
319 }
320 }
321
322 // calculate the well potentials
323 try {
324 updateWellPotentials(reportStepIdx,
325 /*onlyAfterEvent*/true,
326 ebosSimulator_.vanguard().summaryConfig(),
327 local_deferredLogger);
328 } catch ( std::runtime_error& e ) {
329 const std::string msg = "A zero well potential is returned for output purposes. ";
330 local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
331 }
332
333 //update guide rates
334 const auto& comm = ebosSimulator_.vanguard().grid().comm();
335 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
336 std::vector<double> pot(numPhases(), 0.0);
337 const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
338 WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
339 this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
340 std::string exc_msg;
341 auto exc_type = ExceptionType::NONE;
342 // update gpmaint targets
343 if (schedule_[reportStepIdx].has_gpmaint()) {
344 regionalAveragePressureCalculator_->template defineState<ElementContext>(ebosSimulator_);
345 const double dt = ebosSimulator_.timeStepSize();
346 WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup,
347 schedule_, *regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState());
348 }
349 try {
350 // Compute initial well solution for new wells and injectors that change injection type i.e. WAG.
351 for (auto& well : well_container_) {
352 const uint64_t effective_events_mask = ScheduleEvents::WELL_STATUS_CHANGE
353 + ScheduleEvents::INJECTION_TYPE_CHANGED
354 + ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER
355 + ScheduleEvents::NEW_WELL;
356
357 const auto& events = schedule()[reportStepIdx].wellgroup_events();
358 const bool event = report_step_starts_ && events.hasEvent(well->name(), effective_events_mask);
359 const bool dyn_status_change = this->wellState().well(well->name()).status
360 != this->prevWellState().well(well->name()).status;
361
362 if (event || dyn_status_change) {
363 try {
364 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), local_deferredLogger);
365 well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), local_deferredLogger);
366 well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), local_deferredLogger);
367 } catch (const std::exception& e) {
368 const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates";
369 local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
370 }
371 }
372 }
373 }
374 // Catch clauses for all errors setting exc_type and exc_msg
375 OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
376
377 if (exc_type != ExceptionType::NONE) {
378 const std::string msg = "Compute initial well solution for new wells failed. Continue with zero initial rates";
379 local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
380 }
381
382 logAndCheckForExceptionsAndThrow(local_deferredLogger,
383 exc_type, "beginTimeStep() failed: " + exc_msg, terminal_output_, comm);
384
385 }
386
387 template<typename TypeTag>
388 void
389 BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx,
390 const double simulationTime,
391 DeferredLogger& deferred_logger)
392 {
393 const auto& wtest_config = schedule()[timeStepIdx].wtest_config();
394 if (!wtest_config.empty()) { // there is a WTEST request
395 const std::vector<std::string> wellsForTesting = wellTestState()
396 .test_wells(wtest_config, simulationTime);
397 for (const std::string& well_name : wellsForTesting) {
398
399 const Well& wellEcl = schedule().getWell(well_name, timeStepIdx);
400 if (wellEcl.getStatus() == Well::Status::SHUT)
401 continue;
402
403 WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx, deferred_logger);
404 // some preparation before the well can be used
405 well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg_, true);
406
407 double well_efficiency_factor = wellEcl.getEfficiencyFactor();
408 WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx),
409 schedule(), timeStepIdx, well_efficiency_factor);
410
411 well->setWellEfficiencyFactor(well_efficiency_factor);
412 well->setVFPProperties(vfp_properties_.get());
413 well->setGuideRate(&guideRate_);
414
415 well->wellTesting(ebosSimulator_, simulationTime, this->wellState(), this->groupState(), wellTestState(), deferred_logger);
416 }
417 }
418 }
419
420
421
422
423
424 // called at the end of a report step
425 template<typename TypeTag>
426 void
427 BlackoilWellModel<TypeTag>::
428 endReportStep()
429 {
430 // Clear the communication data structures for above values.
431 for (auto&& pinfo : this->local_parallel_well_info_)
432 {
433 pinfo.get().clear();
434 }
435 }
436
437
438
439
440
441 // called at the end of a report step
442 template<typename TypeTag>
443 const SimulatorReportSingle&
444 BlackoilWellModel<TypeTag>::
445 lastReport() const {return last_report_; }
446
447
448
449
450
451 // called at the end of a time step
452 template<typename TypeTag>
453 void
454 BlackoilWellModel<TypeTag>::
455 timeStepSucceeded(const double& simulationTime, const double dt)
456 {
457 this->closed_this_step_.clear();
458
459 // time step is finished and we are not any more at the beginning of an report step
460 report_step_starts_ = false;
461 const int reportStepIdx = ebosSimulator_.episodeIndex();
462
463 DeferredLogger local_deferredLogger;
464 for (const auto& well : well_container_) {
465 if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
466 well->updateWaterThroughput(dt, this->wellState());
467 }
468 }
469 // report well switching
470 for (const auto& well : well_container_) {
471 well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);
472 }
473 // report group switching
474 if (terminal_output_) {
475
476 for (const auto& [name, to] : switched_prod_groups_) {
477 const Group::ProductionCMode& oldControl = this->prevWGState().group_state.production_control(name);
478 std::string from = Group::ProductionCMode2String(oldControl);
479 if (to != from) {
480 std::string msg = " Production Group " + name
481 + " control mode changed from ";
482 msg += from;
483 msg += " to " + to;
484 local_deferredLogger.info(msg);
485 }
486 }
487 for (const auto& [key, to] : switched_inj_groups_) {
488 const std::string& name = key.first;
489 const Opm::Phase& phase = key.second;
490
491 const Group::InjectionCMode& oldControl = this->prevWGState().group_state.injection_control(name, phase);
492 std::string from = Group::InjectionCMode2String(oldControl);
493 if (to != from) {
494 std::string msg = " Injection Group " + name
495 + " control mode changed from ";
496 msg += from;
497 msg += " to " + to;
498 local_deferredLogger.info(msg);
499 }
500 }
501 }
502
503 // update the rate converter with current averages pressures etc in
504 rateConverter_->template defineState<ElementContext>(ebosSimulator_);
505
506 // calculate the well potentials
507 try {
508 updateWellPotentials(reportStepIdx,
509 /*onlyAfterEvent*/false,
510 ebosSimulator_.vanguard().summaryConfig(),
511 local_deferredLogger);
512 } catch ( std::runtime_error& e ) {
513 const std::string msg = "A zero well potential is returned for output purposes. ";
514 local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
515 }
516
517 updateWellTestState(simulationTime, wellTestState());
518
519 // check group sales limits at the end of the timestep
520 const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx);
521 checkGconsaleLimits(fieldGroup, this->wellState(),
522 ebosSimulator_.episodeIndex(), local_deferredLogger);
523
524 this->calculateProductivityIndexValues(local_deferredLogger);
525
526 this->commitWGState();
527
528 const Opm::Parallel::Communication& comm = grid().comm();
529 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
530 if (terminal_output_) {
531 global_deferredLogger.logMessages();
532 }
533
534 //reporting output temperatures
535 this->computeWellTemperature();
536 }
537
538
539 template<typename TypeTag>
540 void
541 BlackoilWellModel<TypeTag>::
542 computeTotalRatesForDof(RateVector& rate,
543 unsigned elemIdx) const
544 {
545 rate = 0;
546
547 if (!is_cell_perforated_[elemIdx])
548 return;
549
550 for (const auto& well : well_container_)
551 well->addCellRates(rate, elemIdx);
552 }
553
554
555 template<typename TypeTag>
556 template <class Context>
557 void
558 BlackoilWellModel<TypeTag>::
559 computeTotalRatesForDof(RateVector& rate,
560 const Context& context,
561 unsigned spaceIdx,
562 unsigned timeIdx) const
563 {
564 rate = 0;
565 int elemIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
566
567 if (!is_cell_perforated_[elemIdx])
568 return;
569
570 for (const auto& well : well_container_)
571 well->addCellRates(rate, elemIdx);
572 }
573
574
575
576 template<typename TypeTag>
577 void
578 BlackoilWellModel<TypeTag>::
579 initializeWellState(const int timeStepIdx,
580 const SummaryState& summaryState)
581 {
582 std::vector<double> cellPressures(this->local_num_cells_, 0.0);
583 ElementContext elemCtx(ebosSimulator_);
584
585 const auto& gridView = ebosSimulator_.vanguard().gridView();
586
587 OPM_BEGIN_PARALLEL_TRY_CATCH();
588 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
589 elemCtx.updatePrimaryStencil(elem);
590 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
591
592 const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState();
593 // copy of get perfpressure in Standard well except for value
594 double& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)];
595 if (Indices::oilEnabled) {
596 perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value();
597 } else if (Indices::waterEnabled) {
598 perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value();
599 } else {
600 perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value();
601 }
602 }
603 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::initializeWellState() failed: ", ebosSimulator_.vanguard().grid().comm());
604
605 this->wellState().init(cellPressures, schedule(), wells_ecl_, local_parallel_well_info_, timeStepIdx,
606 &this->prevWellState(), well_perf_data_,
607 summaryState);
608 }
609
610
611
612
613
614 template<typename TypeTag>
615 void
616 BlackoilWellModel<TypeTag>::
617 createWellContainer(const int time_step)
618 {
619 DeferredLogger local_deferredLogger;
620
621 const int nw = numLocalWells();
622
623 well_container_.clear();
624
625 if (nw > 0) {
626 well_container_.reserve(nw);
627
628 for (int w = 0; w < nw; ++w) {
629 const Well& well_ecl = wells_ecl_[w];
630
631 if (well_ecl.getConnections().empty()) {
632 // No connections in this well. Nothing to do.
633 continue;
634 }
635
636 const std::string& well_name = well_ecl.name();
637 const auto well_status = this->schedule()
638 .getWell(well_name, time_step).getStatus();
639
640 if ((well_ecl.getStatus() == Well::Status::SHUT) ||
641 (well_status == Well::Status::SHUT))
642 {
643 // Due to ACTIONX the well might have been closed behind our back.
644 if (well_ecl.getStatus() != Well::Status::SHUT) {
645 this->closed_this_step_.insert(well_name);
646 this->wellState().shutWell(w);
647 }
648
649 continue;
650 }
651
652 // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
653 if (this->wellTestState().well_is_closed(well_name)) {
654 // TODO: more checking here, to make sure this standard more specific and complete
655 // maybe there is some WCON keywords will not open the well
656 auto& events = this->wellState().well(w).events;
657 if (events.hasEvent(WellState::event_mask)) {
658 if (wellTestState().lastTestTime(well_name) == ebosSimulator_.time()) {
659 // The well was shut this timestep, we are most likely retrying
660 // a timestep without the well in question, after it caused
661 // repeated timestep cuts. It should therefore not be opened,
662 // even if it was new or received new targets this report step.
663 events.clearEvent(WellState::event_mask);
664 } else {
665 wellTestState().open_well(well_name);
666 wellTestState().open_completions(well_name);
667 }
668 }
669 }
670
671 // TODO: should we do this for all kinds of closing reasons?
672 // something like wellTestState().hasWell(well_name)?
673 bool wellIsStopped = false;
674 if (wellTestState().well_is_closed(well_name))
675 {
676 if (well_ecl.getAutomaticShutIn()) {
677 // shut wells are not added to the well container
678 this->wellState().shutWell(w);
679 continue;
680 } else {
681 if (!well_ecl.getAllowCrossFlow()) {
682 // stopped wells where cross flow is not allowed
683 // are not added to the well container
684 this->wellState().shutWell(w);
685 continue;
686 }
687 // stopped wells are added to the container but marked as stopped
688 this->wellState().stopWell(w);
689 wellIsStopped = true;
690 }
691 }
692
693 // If a production well disallows crossflow and its
694 // (prediction type) rate control is zero, then it is effectively shut.
695 if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) {
696 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
697 const auto prod_controls = well_ecl.productionControls(summaryState);
698
699 auto is_zero = [](const double x)
700 {
701 return std::isfinite(x) && !std::isnormal(x);
702 };
703
704 bool zero_rate_control = false;
705 switch (prod_controls.cmode) {
706 case Well::ProducerCMode::ORAT:
707 zero_rate_control = is_zero(prod_controls.oil_rate);
708 break;
709
710 case Well::ProducerCMode::WRAT:
711 zero_rate_control = is_zero(prod_controls.water_rate);
712 break;
713
714 case Well::ProducerCMode::GRAT:
715 zero_rate_control = is_zero(prod_controls.gas_rate);
716 break;
717
718 case Well::ProducerCMode::LRAT:
719 zero_rate_control = is_zero(prod_controls.liquid_rate);
720 break;
721
722 case Well::ProducerCMode::RESV:
723 zero_rate_control = is_zero(prod_controls.resv_rate);
724 break;
725 default:
726 // Might still have zero rate controls, but is pressure controlled.
727 zero_rate_control = false;
728 break;
729 }
730
731 if (zero_rate_control) {
732 // Treat as shut, do not add to container.
733 local_deferredLogger.info(" Well shut due to zero rate control and disallowing crossflow: " + well_ecl.name());
734 this->wellState().shutWell(w);
735 continue;
736 }
737 }
738
739 if (well_status == Well::Status::STOP) {
740 this->wellState().stopWell(w);
741 wellIsStopped = true;
742 }
743
744 well_container_.emplace_back(this->createWellPointer(w, time_step));
745
746 if (wellIsStopped)
747 well_container_.back()->stopWell();
748 }
749 }
750
751 // Collect log messages and print.
752
753 const Opm::Parallel::Communication& comm = grid().comm();
754 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
755 if (terminal_output_) {
756 global_deferredLogger.logMessages();
757 }
758
759 well_container_generic_.clear();
760 for (auto& w : well_container_)
761 well_container_generic_.push_back(w.get());
762 }
763
764
765
766
767
768 template <typename TypeTag>
769 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
770 BlackoilWellModel<TypeTag>::
771 createWellPointer(const int wellID, const int time_step) const
772 {
773 const auto is_multiseg = this->wells_ecl_[wellID].isMultiSegment();
774
775 if (! (this->param_.use_multisegment_well_ && is_multiseg)) {
776 return this->template createTypedWellPointer<StandardWell<TypeTag>>(wellID, time_step);
777 }
778 else {
779 return this->template createTypedWellPointer<MultisegmentWell<TypeTag>>(wellID, time_step);
780 }
781 }
782
783
784
785
786
787 template <typename TypeTag>
788 template <typename WellType>
789 std::unique_ptr<WellType>
790 BlackoilWellModel<TypeTag>::
791 createTypedWellPointer(const int wellID, const int time_step) const
792 {
793 // Use the pvtRegionIdx from the top cell
794 const auto& perf_data = this->well_perf_data_[wellID];
795
796 // Cater for case where local part might have no perforations.
797 const auto pvtreg = perf_data.empty()
798 ? 0 : pvt_region_idx_[perf_data.front().cell_index];
799
800 const auto& parallel_well_info = this->local_parallel_well_info_[wellID].get();
801 const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
802
803 return std::make_unique<WellType>(this->wells_ecl_[wellID],
804 parallel_well_info,
805 time_step,
806 this->param_,
807 *this->rateConverter_,
808 global_pvtreg,
809 this->numComponents(),
810 this->numPhases(),
811 wellID,
812 perf_data);
813 }
814
815
816
817
818
819 template<typename TypeTag>
820 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
821 BlackoilWellModel<TypeTag>::
822 createWellForWellTest(const std::string& well_name,
823 const int report_step,
824 DeferredLogger& deferred_logger) const
825 {
826 // Finding the location of the well in wells_ecl
827 const int nw_wells_ecl = wells_ecl_.size();
828 int index_well_ecl = 0;
829 for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
830 if (well_name == wells_ecl_[index_well_ecl].name()) {
831 break;
832 }
833 }
834 // It should be able to find in wells_ecl.
835 if (index_well_ecl == nw_wells_ecl) {
836 OPM_DEFLOG_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ", deferred_logger);
837 }
838
839 return this->createWellPointer(index_well_ecl, report_step);
840 }
841
842
843
844
845
846 template<typename TypeTag>
847 void
848 BlackoilWellModel<TypeTag>::
849 assemble(const int iterationIdx,
850 const double dt)
851 {
852
853 DeferredLogger local_deferredLogger;
854 if (this->glift_debug) {
855 const std::string msg = fmt::format(
856 "assemble() : iteration {}" , iterationIdx);
857 gliftDebug(msg, local_deferredLogger);
858 }
859 last_report_ = SimulatorReportSingle();
860 Dune::Timer perfTimer;
861 perfTimer.start();
862
863 if ( ! wellsActive() ) {
864 return;
865 }
866
867 updatePerforationIntensiveQuantities();
868
869 if (iterationIdx == 0) {
870 // try-catch is needed here as updateWellControls
871 // contains global communication and has either to
872 // be reached by all processes or all need to abort
873 // before.
874 OPM_BEGIN_PARALLEL_TRY_CATCH();
875 {
876 calculateExplicitQuantities(local_deferredLogger);
877 prepareTimeStep(local_deferredLogger);
878 }
879 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed (It=0): ",
880 terminal_output_, grid().comm());
881 }
882
883 const bool well_group_control_changed = assembleImpl(iterationIdx, dt, 0, local_deferredLogger);
884
885 // if group or well control changes we don't consider the
886 // case converged
887 last_report_.well_group_control_changed = well_group_control_changed;
888 last_report_.assemble_time_well += perfTimer.stop();
889 }
890
891
892
893
894
895 template<typename TypeTag>
896 bool
897 BlackoilWellModel<TypeTag>::
898 assembleImpl(const int iterationIdx,
899 const double dt,
900 const std::size_t recursion_level,
901 DeferredLogger& local_deferredLogger)
902 {
903
904 auto [well_group_control_changed, network_changed, network_imbalance] = updateWellControls(local_deferredLogger);
905
906 bool alq_updated = false;
907 OPM_BEGIN_PARALLEL_TRY_CATCH();
908 {
909 // Set the well primary variables based on the value of well solutions
910 initPrimaryVariablesEvaluation();
911
912 alq_updated = maybeDoGasLiftOptimize(local_deferredLogger);
913 assembleWellEq(dt, local_deferredLogger);
914 }
915 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed: ",
916 terminal_output_, grid().comm());
917
918 //update guide rates
919 const int reportStepIdx = ebosSimulator_.episodeIndex();
920 if (alq_updated || guideRateUpdateIsNeeded(reportStepIdx)) {
921 const double simulationTime = ebosSimulator_.time();
922 const auto& comm = ebosSimulator_.vanguard().grid().comm();
923 const auto& summaryState = ebosSimulator_.vanguard().summaryState();
924 std::vector<double> pot(numPhases(), 0.0);
925 const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
926 WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
927 this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
928 }
929
930 // Maybe do a recursive call to iterate network and well controls.
931 if (network_changed) {
932 if (shouldBalanceNetwork(reportStepIdx, iterationIdx)) {
933 const auto& balance = schedule()[reportStepIdx].network_balance();
934 // Iterate if not converged, and number of iterations is not yet max (NETBALAN item 3).
935 if (recursion_level < balance.pressure_max_iter() && network_imbalance > balance.pressure_tolerance()) {
936 well_group_control_changed = assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
937 }
938 }
939 }
940 return well_group_control_changed;
941 }
942
943
944
945
946
947 template<typename TypeTag>
948 bool
949 BlackoilWellModel<TypeTag>::
950 maybeDoGasLiftOptimize(DeferredLogger& deferred_logger)
951 {
952 bool do_glift_optimization = false;
953 int num_wells_changed = 0;
954 const double simulation_time = ebosSimulator_.time();
955 const double min_wait = ebosSimulator_.vanguard().schedule().glo(ebosSimulator_.episodeIndex()).min_wait();
956 // We only optimize if a min_wait time has past.
957 // If all_newton is true we still want to optimize several times pr timestep
958 // i.e. we also optimize if check simulation_time == last_glift_opt_time_
959 // that is when the last_glift_opt_time is already updated with the current time step
960 if ( simulation_time == last_glift_opt_time_ || simulation_time >= (last_glift_opt_time_ + min_wait)) {
961 do_glift_optimization = true;
962 last_glift_opt_time_ = simulation_time;
963 }
964
965 if (do_glift_optimization) {
966 GLiftOptWells glift_wells;
967 GLiftProdWells prod_wells;
968 GLiftWellStateMap state_map;
969 // NOTE: To make GasLiftGroupInfo (see below) independent of the TypeTag
970 // associated with *this (i.e. BlackoilWellModel<TypeTag>) we observe
971 // that GasLiftGroupInfo's only dependence on *this is that it needs to
972 // access the eclipse Wells in the well container (the eclipse Wells
973 // themselves are independent of the TypeTag).
974 // Hence, we extract them from the well container such that we can pass
975 // them to the GasLiftGroupInfo constructor.
976 GLiftEclWells ecl_well_map;
977 initGliftEclWellMap(ecl_well_map);
978 GasLiftGroupInfo group_info {
979 ecl_well_map,
980 ebosSimulator_.vanguard().schedule(),
981 ebosSimulator_.vanguard().summaryState(),
982 ebosSimulator_.episodeIndex(),
983 ebosSimulator_.model().newtonMethod().numIterations(),
984 phase_usage_,
985 deferred_logger,
986 this->wellState(),
987 ebosSimulator_.vanguard().grid().comm(),
988 this->glift_debug
989 };
990 group_info.initialize();
991 gasLiftOptimizationStage1(
992 deferred_logger, prod_wells, glift_wells, group_info, state_map);
993 gasLiftOptimizationStage2(
994 deferred_logger, prod_wells, glift_wells, state_map,
995 ebosSimulator_.episodeIndex());
996 if (this->glift_debug) gliftDebugShowALQ(deferred_logger);
997 num_wells_changed = glift_wells.size();
998 }
999 num_wells_changed = this->comm_.sum(num_wells_changed);
1000 return num_wells_changed > 0;
1001 }
1002
1003 template<typename TypeTag>
1004 void
1005 BlackoilWellModel<TypeTag>::
1006 gasLiftOptimizationStage1(DeferredLogger& deferred_logger,
1007 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
1008 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map)
1009 {
1010 auto comm = ebosSimulator_.vanguard().grid().comm();
1011 int num_procs = comm.size();
1012 // NOTE: Gas lift optimization stage 1 seems to be difficult
1013 // to do in parallel since the wells are optimized on different
1014 // processes and each process needs to know the current ALQ allocated
1015 // to each group it is a memeber of in order to check group limits and avoid
1016 // allocating more ALQ than necessary. (Surplus ALQ is removed in
1017 // stage 2). In stage1, as each well is adding ALQ, the current group ALQ needs
1018 // to be communicated to the other processes. But there is no common
1019 // synchronization point that all process will reach in the
1020 // runOptimizeLoop_() in GasLiftSingleWell.cpp.
1021 //
1022 // TODO: Maybe a better solution could be invented by distributing
1023 // wells according to certain parent groups. Then updated group rates
1024 // might not have to be communicated to the other processors.
1025
1026 // Currently, the best option seems to be to run this part sequentially
1027 // (not in parallel).
1028 //
1029 // TODO: The simplest approach seems to be if a) one process could take
1030 // ownership of all the wells (the union of all the wells in the
1031 // well_container_ of each process) then this process could do the
1032 // optimization, while the other processes could wait for it to
1033 // finish (e.g. comm.barrier()), or alternatively, b) if all
1034 // processes could take ownership of all the wells. Then there
1035 // would be no need for synchronization here..
1036 //
1037 for (int i = 0; i< num_procs; i++) {
1038 int num_rates_to_sync = 0; // communication variable
1039 GLiftSyncGroups groups_to_sync;
1040 if (comm.rank() == i) {
1041 // Run stage1: Optimize single wells while also checking group limits
1042 for (const auto& well : well_container_) {
1043 // NOTE: Only the wells in "group_info" needs to be optimized
1044 if (group_info.hasWell(well->name())) {
1045 gasLiftOptimizationStage1SingleWell(
1046 well.get(), deferred_logger, prod_wells, glift_wells,
1047 group_info, state_map, groups_to_sync
1048 );
1049 }
1050 }
1051 num_rates_to_sync = groups_to_sync.size();
1052 }
1053 // Since "group_info" is not used in stage2, there is no need to
1054 // communicate rates if this is the last iteration...
1055 if (i == (num_procs - 1))
1056 break;
1057 num_rates_to_sync = comm.sum(num_rates_to_sync);
1058 if (num_rates_to_sync > 0) {
1059 std::vector<int> group_indexes;
1060 group_indexes.reserve(num_rates_to_sync);
1061 std::vector<double> group_alq_rates;
1062 group_alq_rates.reserve(num_rates_to_sync);
1063 std::vector<double> group_oil_rates;
1064 group_oil_rates.reserve(num_rates_to_sync);
1065 std::vector<double> group_gas_rates;
1066 group_gas_rates.reserve(num_rates_to_sync);
1067 std::vector<double> group_water_rates;
1068 group_water_rates.reserve(num_rates_to_sync);
1069 if (comm.rank() == i) {
1070 for (auto idx : groups_to_sync) {
1071 auto [oil_rate, gas_rate, water_rate, alq] = group_info.getRates(idx);
1072 group_indexes.push_back(idx);
1073 group_oil_rates.push_back(oil_rate);
1074 group_gas_rates.push_back(gas_rate);
1075 group_water_rates.push_back(water_rate);
1076 group_alq_rates.push_back(alq);
1077 }
1078 } else {
1079 group_indexes.resize(num_rates_to_sync);
1080 group_oil_rates.resize(num_rates_to_sync);
1081 group_gas_rates.resize(num_rates_to_sync);
1082 group_water_rates.resize(num_rates_to_sync);
1083 group_alq_rates.resize(num_rates_to_sync);
1084 }
1085 // TODO: We only need to broadcast to processors with index
1086 // j > i since we do not use the "group_info" in stage 2. In
1087 // this case we should use comm.send() and comm.receive()
1088 // instead of comm.broadcast() to communicate with specific
1089 // processes, and these processes only need to receive the
1090 // data if they are going to check the group rates in stage1
1091 // Another similar idea is to only communicate the rates to
1092 // process j = i + 1
1093#if HAVE_MPI
1094 EclMpiSerializer ser(comm);
1095 ser.broadcast(i, group_indexes, group_oil_rates,
1096 group_gas_rates, group_water_rates, group_alq_rates);
1097#endif
1098 if (comm.rank() != i) {
1099 for (int j=0; j<num_rates_to_sync; j++) {
1100 group_info.updateRate(group_indexes[j],
1101 group_oil_rates[j], group_gas_rates[j], group_water_rates[j], group_alq_rates[j]);
1102 }
1103 }
1104 if (this->glift_debug) {
1105 int counter = 0;
1106 if (comm.rank() == i) {
1107 counter = this->wellState().gliftGetDebugCounter();
1108 }
1109 counter = comm.sum(counter);
1110 if (comm.rank() != i) {
1111 this->wellState().gliftSetDebugCounter(counter);
1112 }
1113 }
1114 }
1115 }
1116 }
1117
1118 // NOTE: this method cannot be const since it passes this->wellState()
1119 // (see below) to the GasLiftSingleWell constructor which accepts WellState
1120 // as a non-const reference..
1121 template<typename TypeTag>
1122 void
1123 BlackoilWellModel<TypeTag>::
1124 gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
1125 DeferredLogger& deferred_logger,
1126 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
1127 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
1128 GLiftSyncGroups& sync_groups)
1129 {
1130 const auto& summary_state = ebosSimulator_.vanguard().summaryState();
1131 std::unique_ptr<GasLiftSingleWell> glift
1132 = std::make_unique<GasLiftSingleWell>(
1133 *well, ebosSimulator_, summary_state,
1134 deferred_logger, this->wellState(), this->groupState(),
1135 group_info, sync_groups, this->comm_, this->glift_debug);
1136 auto state = glift->runOptimize(
1137 ebosSimulator_.model().newtonMethod().numIterations());
1138 if (state) {
1139 state_map.insert({well->name(), std::move(state)});
1140 glift_wells.insert({well->name(), std::move(glift)});
1141 return;
1142 }
1143 prod_wells.insert({well->name(), well});
1144 }
1145
1146
1147 template<typename TypeTag>
1148 void
1149 BlackoilWellModel<TypeTag>::
1150 initGliftEclWellMap(GLiftEclWells &ecl_well_map)
1151 {
1152 for ( const auto& well: well_container_ ) {
1153 ecl_well_map.try_emplace(
1154 well->name(), &(well->wellEcl()), well->indexOfWell());
1155 }
1156 }
1157
1158
1159 template<typename TypeTag>
1160 void
1161 BlackoilWellModel<TypeTag>::
1162 assembleWellEq(const double dt, DeferredLogger& deferred_logger)
1163 {
1164 for (auto& well : well_container_) {
1165 well->assembleWellEq(ebosSimulator_, dt, this->wellState(), this->groupState(), deferred_logger);
1166 }
1167 }
1168
1169 template<typename TypeTag>
1170 void
1171 BlackoilWellModel<TypeTag>::
1172 apply( BVector& r) const
1173 {
1174 for (auto& well : well_container_) {
1175 well->apply(r);
1176 }
1177 }
1178
1179
1180 // Ax = A x - C D^-1 B x
1181 template<typename TypeTag>
1182 void
1183 BlackoilWellModel<TypeTag>::
1184 apply(const BVector& x, BVector& Ax) const
1185 {
1186 for (auto& well : well_container_) {
1187 well->apply(x, Ax);
1188 }
1189 }
1190
1191 template<typename TypeTag>
1192 void
1193 BlackoilWellModel<TypeTag>::
1194 getWellContributions(WellContributions& wellContribs) const
1195 {
1196 // prepare for StandardWells
1197 wellContribs.setBlockSize(StandardWell<TypeTag>::Indices::numEq, StandardWell<TypeTag>::numStaticWellEq);
1198
1199 for(unsigned int i = 0; i < well_container_.size(); i++){
1200 auto& well = well_container_[i];
1201 std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1202 if (derived) {
1203 unsigned int numBlocks;
1204 derived->getNumBlocks(numBlocks);
1205 wellContribs.addNumBlocks(numBlocks);
1206 }
1207 }
1208
1209 // allocate memory for data from StandardWells
1210 wellContribs.alloc();
1211
1212 for(unsigned int i = 0; i < well_container_.size(); i++){
1213 auto& well = well_container_[i];
1214 // maybe WellInterface could implement addWellContribution()
1215 auto derived_std = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);
1216 if (derived_std) {
1217 derived_std->addWellContribution(wellContribs);
1218 } else {
1219 auto derived_ms = std::dynamic_pointer_cast<MultisegmentWell<TypeTag> >(well);
1220 if (derived_ms) {
1221 derived_ms->addWellContribution(wellContribs);
1222 } else {
1223 OpmLog::warning("Warning unknown type of well");
1224 }
1225 }
1226 }
1227 }
1228
1229 // Ax = Ax - alpha * C D^-1 B x
1230 template<typename TypeTag>
1231 void
1232 BlackoilWellModel<TypeTag>::
1233 applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
1234 {
1235 if (this->well_container_.empty()) {
1236 return;
1237 }
1238
1239 if( scaleAddRes_.size() != Ax.size() ) {
1240 scaleAddRes_.resize( Ax.size() );
1241 }
1242
1243 scaleAddRes_ = 0.0;
1244 // scaleAddRes_ = - C D^-1 B x
1245 apply( x, scaleAddRes_ );
1246 // Ax = Ax + alpha * scaleAddRes_
1247 Ax.axpy( alpha, scaleAddRes_ );
1248 }
1249
1250 template<typename TypeTag>
1251 void
1252 BlackoilWellModel<TypeTag>::
1253 addWellContributions(SparseMatrixAdapter& jacobian) const
1254 {
1255 for ( const auto& well: well_container_ ) {
1256 well->addWellContributions(jacobian);
1257 }
1258 }
1259
1260 template<typename TypeTag>
1261 void
1262 BlackoilWellModel<TypeTag>::
1263 addWellPressureEquations(PressureMatrix& jacobian, const BVector& weights,const bool use_well_weights) const
1264 {
1265 int nw = this->numLocalWellsEnd();
1266 int rdofs = local_num_cells_;
1267 for ( int i = 0; i < nw; i++ ){
1268 int wdof = rdofs + i;
1269 jacobian[wdof][wdof] = 1.0;// better scaling ?
1270 }
1271
1272 for ( const auto& well : well_container_ ) {
1273 well->addWellPressureEquations(jacobian, weights, pressureVarIndex, use_well_weights, this->wellState());
1274 }
1275 }
1276
1277 template<typename TypeTag>
1278 int
1279 BlackoilWellModel<TypeTag>::
1280 numLocalWellsEnd() const
1281 {
1282 auto w = schedule().getWellsatEnd();
1283 w.erase(std::remove_if(w.begin(), w.end(), not_on_process_), w.end());
1284 return w.size();
1285 }
1286
1287 template<typename TypeTag>
1288 std::vector<std::vector<int>>
1289 BlackoilWellModel<TypeTag>::
1290 getMaxWellConnections() const
1291 {
1292 std::vector<std::vector<int>> wells;
1293
1294 auto schedule_wells = schedule().getWellsatEnd();
1295 schedule_wells.erase(std::remove_if(schedule_wells.begin(), schedule_wells.end(), not_on_process_), schedule_wells.end());
1296 wells.reserve(schedule_wells.size());
1297
1298 // initialize the additional cell connections introduced by wells.
1299 for ( const auto& well : schedule_wells )
1300 {
1301 std::vector<int> compressed_well_perforations;
1302 // All possible completions of the well
1303 const auto& completionSet = well.getConnections();
1304 compressed_well_perforations.reserve(completionSet.size());
1305
1306 for (const auto& connection: well.getConnections())
1307 {
1308 const int compressed_idx = compressedIndexForInterior(connection.global_index());
1309 if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells.
1310 {
1311 compressed_well_perforations.push_back(compressed_idx);
1312 }
1313 }
1314
1315 // also include wells with no perforations in case
1316 std::sort(compressed_well_perforations.begin(),
1317 compressed_well_perforations.end());
1318
1319 wells.push_back(compressed_well_perforations);
1320 }
1321 return wells;
1322 }
1323
1324 template<typename TypeTag>
1325 void
1326 BlackoilWellModel<TypeTag>::
1327 addWellPressureEquationsStruct(PressureMatrix& jacobian) const
1328 {
1329 int nw = this->numLocalWellsEnd();
1330 int rdofs = local_num_cells_;
1331 for(int i=0; i < nw; i++){
1332 int wdof = rdofs + i;
1333 jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
1334 }
1335 std::vector<std::vector<int>> wellconnections = getMaxWellConnections();
1336 for(int i=0; i < nw; i++){
1337 const auto& perfcells = wellconnections[i];
1338 for(int perfcell : perfcells){
1339 int wdof = rdofs + i;
1340 jacobian.entry(wdof,perfcell) = 0.0;
1341 jacobian.entry(perfcell, wdof) = 0.0;
1342 }
1343 }
1344 }
1345
1346 template<typename TypeTag>
1347 int
1348 BlackoilWellModel<TypeTag>::
1349 numLocalNonshutWells() const
1350 {
1351 return well_container_.size();
1352 }
1353
1354 template<typename TypeTag>
1355 void
1356 BlackoilWellModel<TypeTag>::
1357 recoverWellSolutionAndUpdateWellState(const BVector& x)
1358 {
1359
1360 DeferredLogger local_deferredLogger;
1361 OPM_BEGIN_PARALLEL_TRY_CATCH();
1362 {
1363 for (auto& well : well_container_) {
1364 well->recoverWellSolutionAndUpdateWellState(x, this->wellState(), local_deferredLogger);
1365 }
1366
1367 }
1368 OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger,
1369 "recoverWellSolutionAndUpdateWellState() failed: ",
1370 terminal_output_, ebosSimulator_.vanguard().grid().comm());
1371
1372 }
1373
1374
1375
1376
1377 template<typename TypeTag>
1378 void
1379 BlackoilWellModel<TypeTag>::
1380 initPrimaryVariablesEvaluation() const
1381 {
1382 for (auto& well : well_container_) {
1383 well->initPrimaryVariablesEvaluation();
1384 }
1385 }
1386
1387
1388
1389
1390
1391
1392 template<typename TypeTag>
1393 ConvergenceReport
1394 BlackoilWellModel<TypeTag>::
1395 getWellConvergence(const std::vector<Scalar>& B_avg, bool checkWellGroupControls) const
1396 {
1397
1398 DeferredLogger local_deferredLogger;
1399 // Get global (from all processes) convergence report.
1400 ConvergenceReport local_report;
1401 const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1402 for (const auto& well : well_container_) {
1403 if (well->isOperableAndSolvable() || well->wellIsStopped()) {
1404 local_report += well->getWellConvergence(this->wellState(), B_avg, local_deferredLogger, iterationIdx > param_.strict_outer_iter_wells_ );
1405 } else {
1406 ConvergenceReport report;
1407 using CR = ConvergenceReport;
1408 report.setWellFailed({CR::WellFailure::Type::Unsolvable, CR::Severity::Normal, -1, well->name()});
1409 local_report += report;
1410 }
1411 }
1412
1413 const Opm::Parallel::Communication comm = grid().comm();
1414 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
1415 ConvergenceReport report = gatherConvergenceReport(local_report, comm);
1416
1417 // the well_group_control_changed info is already communicated
1418 if (checkWellGroupControls) {
1419 report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed);
1420 }
1421
1422 if (terminal_output_) {
1423 global_deferredLogger.logMessages();
1424 }
1425 // Log debug messages for NaN or too large residuals.
1426 if (terminal_output_) {
1427 for (const auto& f : report.wellFailures()) {
1428 if (f.severity() == ConvergenceReport::Severity::NotANumber) {
1429 OpmLog::debug("NaN residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1430 } else if (f.severity() == ConvergenceReport::Severity::TooLarge) {
1431 OpmLog::debug("Too large residual found with phase " + std::to_string(f.phase()) + " for well " + f.wellName());
1432 }
1433 }
1434 }
1435 return report;
1436 }
1437
1438
1439
1440
1441
1442 template<typename TypeTag>
1443 void
1445 calculateExplicitQuantities(DeferredLogger& deferred_logger) const
1446 {
1447 // TODO: checking isOperableAndSolvable() ?
1448 for (auto& well : well_container_) {
1449 well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), deferred_logger);
1450 }
1451 }
1452
1453
1454
1455
1456
1457 template<typename TypeTag>
1458 bool
1460 shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const
1461 {
1462 const auto& balance = schedule()[reportStepIdx].network_balance();
1463 if (balance.mode() == Network::Balance::CalcMode::TimeStepStart) {
1464 return iterationIdx == 0;
1465 } else if (balance.mode() == Network::Balance::CalcMode::NUPCOL) {
1466 const int nupcol = schedule()[reportStepIdx].nupcol();
1467 return iterationIdx < nupcol;
1468 } else {
1469 // We do not support any other rebalancing modes,
1470 // i.e. TimeInterval based rebalancing is not available.
1471 // This should be warned about elsewhere, so we choose to
1472 // avoid spamming with a warning here.
1473 return false;
1474 }
1475 }
1476
1477
1478
1479
1480
1481 template<typename TypeTag>
1482 std::tuple<bool, bool, double>
1483 BlackoilWellModel<TypeTag>::
1484 updateWellControls(DeferredLogger& deferred_logger)
1485 {
1486 // Even if there are no wells active locally, we cannot
1487 // return as the DeferredLogger uses global communication.
1488 // For no well active globally we simply return.
1489 if( !wellsActive() ) return { false, false, 0.0 };
1490
1491 const int episodeIdx = ebosSimulator_.episodeIndex();
1492 const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
1493 const auto& comm = ebosSimulator_.vanguard().grid().comm();
1494 updateAndCommunicateGroupData(episodeIdx, iterationIdx);
1495
1496 const auto [local_network_changed, local_network_imbalance]
1497 = shouldBalanceNetwork(episodeIdx, iterationIdx) ?
1498 updateNetworkPressures(episodeIdx) : std::make_pair(false, 0.0);
1499 const bool network_changed = comm.sum(local_network_changed);
1500 const double network_imbalance = comm.max(local_network_imbalance);
1501
1502 bool changed_well_group = false;
1503 // Check group individual constraints.
1504 const int nupcol = schedule()[episodeIdx].nupcol();
1505 // don't switch group control when iterationIdx > nupcol
1506 // to avoid oscilations between group controls
1507 if (iterationIdx <= nupcol) {
1508 const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
1509 changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
1510 }
1511 // Check wells' group constraints and communicate.
1512 bool changed_well_to_group = false;
1513 for (const auto& well : well_container_) {
1514 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
1515 const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1516 if (changed_well) {
1517 changed_well_to_group = changed_well || changed_well_to_group;
1518 }
1519 }
1520
1521 changed_well_to_group = comm.sum(changed_well_to_group);
1522 if (changed_well_to_group) {
1523 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1524 changed_well_group = true;
1525 }
1526
1527 // Check individual well constraints and communicate.
1528 bool changed_well_individual = false;
1529 for (const auto& well : well_container_) {
1530 const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
1531 const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
1532 if (changed_well) {
1533 changed_well_individual = changed_well || changed_well_individual;
1534 }
1535 }
1536 changed_well_individual = comm.sum(changed_well_individual);
1537 if (changed_well_individual) {
1538 updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
1539 changed_well_group = true;
1540 }
1541
1542 // update wsolvent fraction for REIN wells
1543 const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
1544 updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState());
1545
1546 return { changed_well_group, network_changed, network_imbalance };
1547 }
1548
1549
1550 template<typename TypeTag>
1551 void
1552 BlackoilWellModel<TypeTag>::
1553 updateAndCommunicate(const int reportStepIdx,
1554 const int iterationIdx,
1555 DeferredLogger& deferred_logger)
1556 {
1557 updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
1558 // if a well or group change control it affects all wells that are under the same group
1559 for (const auto& well : well_container_) {
1560 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
1561 }
1562 updateAndCommunicateGroupData(reportStepIdx, iterationIdx);
1563 }
1564
1565 template<typename TypeTag>
1566 bool
1567 BlackoilWellModel<TypeTag>::
1568 updateGroupControls(const Group& group,
1569 DeferredLogger& deferred_logger,
1570 const int reportStepIdx,
1571 const int iterationIdx)
1572 {
1573 bool changed = false;
1574 bool changed_hc = checkGroupHigherConstraints( group, deferred_logger, reportStepIdx);
1575 if (changed_hc) {
1576 changed = true;
1577 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1578 }
1579 bool changed_individual = updateGroupIndividualControl( group, deferred_logger, reportStepIdx);
1580 if (changed_individual) {
1581 changed = true;
1582 updateAndCommunicate(reportStepIdx, iterationIdx, deferred_logger);
1583 }
1584 // call recursively down the group hierarchy
1585 for (const std::string& groupName : group.groups()) {
1586 bool changed_this = updateGroupControls( schedule().getGroup(groupName, reportStepIdx), deferred_logger, reportStepIdx,iterationIdx);
1587 changed = changed || changed_this;
1588 }
1589 return changed;
1590 }
1591
1592 template<typename TypeTag>
1593 void
1595 updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
1596 {
1597 DeferredLogger local_deferredLogger;
1598 for (const auto& well : well_container_) {
1599 const auto& wname = well->name();
1600 const auto wasClosed = wellTestState.well_is_closed(wname);
1601 well->checkWellOperability(ebosSimulator_, this->wellState(), local_deferredLogger);
1602 well->updateWellTestState(this->wellState().well(wname), simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
1603
1604 if (!wasClosed && wellTestState.well_is_closed(wname)) {
1605 this->closed_this_step_.insert(wname);
1606 }
1607 }
1608
1609 const Opm::Parallel::Communication comm = grid().comm();
1610 DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
1611 if (terminal_output_) {
1612 global_deferredLogger.logMessages();
1613 }
1614 }
1615
1616
1617 template<typename TypeTag>
1618 void
1619 BlackoilWellModel<TypeTag>::computePotentials(const std::size_t widx,
1620 const WellState& well_state_copy,
1621 std::string& exc_msg,
1622 ExceptionType::ExcEnum& exc_type,
1623 DeferredLogger& deferred_logger)
1624 {
1625 const int np = numPhases();
1626 std::vector<double> potentials;
1627 const auto& well= well_container_[widx];
1628 try {
1629 well->computeWellPotentials(ebosSimulator_, well_state_copy, potentials, deferred_logger);
1630 }
1631 // catch all possible exception and store type and message.
1632 OPM_PARALLEL_CATCH_CLAUSE(exc_type, exc_msg);
1633 // Store it in the well state
1634 // potentials is resized and set to zero in the beginning of well->ComputeWellPotentials
1635 // and updated only if sucessfull. i.e. the potentials are zero for exceptions
1636 auto& ws = this->wellState().well(well->indexOfWell());
1637 for (int p = 0; p < np; ++p) {
1638 // make sure the potentials are positive
1639 ws.well_potentials[p] = std::max(0.0, potentials[p]);
1640 }
1641 }
1642
1643
1644
1645 template <typename TypeTag>
1646 void
1647 BlackoilWellModel<TypeTag>::
1648 calculateProductivityIndexValues(DeferredLogger& deferred_logger)
1649 {
1650 for (const auto& wellPtr : this->well_container_) {
1651 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
1652 }
1653 }
1654
1655
1656
1657
1658
1659 template <typename TypeTag>
1660 void
1661 BlackoilWellModel<TypeTag>::
1662 calculateProductivityIndexValuesShutWells(const int reportStepIdx,
1663 DeferredLogger& deferred_logger)
1664 {
1665 // For the purpose of computing PI/II values, it is sufficient to
1666 // construct StandardWell instances only. We don't need to form
1667 // well objects that honour the 'isMultisegment()' flag of the
1668 // corresponding "this->wells_ecl_[shutWell]".
1669
1670 for (const auto& shutWell : this->local_shut_wells_) {
1671 if (this->wells_ecl_[shutWell].getConnections().empty()) {
1672 // No connections in this well. Nothing to do.
1673 continue;
1674 }
1675
1676 auto wellPtr = this->template createTypedWellPointer
1677 <StandardWell<TypeTag>>(shutWell, reportStepIdx);
1678
1679 wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
1680 this->local_num_cells_, this->B_avg_, true);
1681
1682 this->calculateProductivityIndexValues(wellPtr.get(), deferred_logger);
1683 }
1684 }
1685
1686
1687
1688
1689
1690 template <typename TypeTag>
1691 void
1692 BlackoilWellModel<TypeTag>::
1693 calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
1694 DeferredLogger& deferred_logger)
1695 {
1696 wellPtr->updateProductivityIndex(this->ebosSimulator_,
1697 this->prod_index_calc_[wellPtr->indexOfWell()],
1698 this->wellState(),
1699 deferred_logger);
1700 }
1701
1702
1703
1704
1705
1706 template<typename TypeTag>
1707 void
1708 BlackoilWellModel<TypeTag>::
1709 prepareTimeStep(DeferredLogger& deferred_logger)
1710 {
1711 for (const auto& well : well_container_) {
1712 auto& events = this->wellState().well(well->indexOfWell()).events;
1713 if (events.hasEvent(WellState::event_mask)) {
1714 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
1715 well->updatePrimaryVariables(this->wellState(), deferred_logger);
1716 well->initPrimaryVariablesEvaluation();
1717 // There is no new well control change input within a report step,
1718 // so next time step, the well does not consider to have effective events anymore.
1719 events.clearEvent(WellState::event_mask);
1720 }
1721 // solve the well equation initially to improve the initial solution of the well model
1722 if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
1723 try {
1724 well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), deferred_logger);
1725 } catch (const std::exception& e) {
1726 const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the privious rates";
1727 deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
1728 }
1729 }
1730 }
1731 updatePrimaryVariables(deferred_logger);
1732 }
1733
1734 template<typename TypeTag>
1735 void
1736 BlackoilWellModel<TypeTag>::
1737 updateAverageFormationFactor()
1738 {
1739 std::vector< Scalar > B_avg(numComponents(), Scalar() );
1740 const auto& grid = ebosSimulator_.vanguard().grid();
1741 const auto& gridView = grid.leafGridView();
1742 ElementContext elemCtx(ebosSimulator_);
1743
1744 OPM_BEGIN_PARALLEL_TRY_CATCH();
1745 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1746 elemCtx.updatePrimaryStencil(elem);
1747 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1748
1749 const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
1750 const auto& fs = intQuants.fluidState();
1751
1752 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1753 {
1754 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1755 continue;
1756 }
1757
1758 const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1759 auto& B = B_avg[ compIdx ];
1760
1761 B += 1 / fs.invB(phaseIdx).value();
1762 }
1763 if constexpr (has_solvent_) {
1764 auto& B = B_avg[solventSaturationIdx];
1765 B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
1766 }
1767 }
1768 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updateAverageFormationFactor() failed: ", grid.comm())
1769
1770 // compute global average
1771 grid.comm().sum(B_avg.data(), B_avg.size());
1772 for (auto& bval : B_avg)
1773 {
1774 bval /= global_num_cells_;
1775 }
1776 B_avg_ = B_avg;
1777 }
1778
1779
1780
1781
1782
1783 template<typename TypeTag>
1784 void
1785 BlackoilWellModel<TypeTag>::
1786 updatePrimaryVariables(DeferredLogger& deferred_logger)
1787 {
1788 for (const auto& well : well_container_) {
1789 well->updatePrimaryVariables(this->wellState(), deferred_logger);
1790 }
1791 }
1792
1793 template<typename TypeTag>
1794 void
1795 BlackoilWellModel<TypeTag>::extractLegacyCellPvtRegionIndex_()
1796 {
1797 const auto& grid = ebosSimulator_.vanguard().grid();
1798 const auto& eclProblem = ebosSimulator_.problem();
1799 const unsigned numCells = grid.size(/*codim=*/0);
1800
1801 pvt_region_idx_.resize(numCells);
1802 for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
1803 pvt_region_idx_[cellIdx] =
1804 eclProblem.pvtRegionIndex(cellIdx);
1805 }
1806 }
1807
1808 // The number of components in the model.
1809 template<typename TypeTag>
1810 int
1811 BlackoilWellModel<TypeTag>::numComponents() const
1812 {
1813 // The numComponents here does not reflect the actual number of the components in the system.
1814 // It more or less reflects the number of mass conservation equations for the well equations.
1815 // For example, in the current formulation, we do not have the polymer conservation equation
1816 // in the well equations. As a result, for an oil-water-polymer system, this function will return 2.
1817 // In some way, it makes this function appear to be confusing from its name, and we need
1818 // to revisit/revise this function again when extending the variants of system that flow can simulate.
1819 if (numPhases() < 3) {
1820 return numPhases();
1821 }
1822 int numComp = FluidSystem::numComponents;
1823 if constexpr (has_solvent_) {
1824 numComp ++;
1825 }
1826
1827 return numComp;
1828 }
1829
1830 template<typename TypeTag>
1831 void
1832 BlackoilWellModel<TypeTag>::extractLegacyDepth_()
1833 {
1834 const auto& eclProblem = ebosSimulator_.problem();
1835 depth_.resize(local_num_cells_);
1836 for (unsigned cellIdx = 0; cellIdx < local_num_cells_; ++cellIdx) {
1837 depth_[cellIdx] = eclProblem.dofCenterDepth(cellIdx);
1838 }
1839 }
1840
1841 template<typename TypeTag>
1842 void
1843 BlackoilWellModel<TypeTag>::
1844 updatePerforationIntensiveQuantities()
1845 {
1846 ElementContext elemCtx(ebosSimulator_);
1847 const auto& gridView = ebosSimulator_.gridView();
1848
1849 OPM_BEGIN_PARALLEL_TRY_CATCH();
1850 for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
1851 elemCtx.updatePrimaryStencil(elem);
1852 int elemIdx = elemCtx.globalSpaceIndex(0, 0);
1853
1854 if (!is_cell_perforated_[elemIdx]) {
1855 continue;
1856 }
1857 elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
1858 }
1859 OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updatePerforationIntensiveQuantities() failed: ", ebosSimulator_.vanguard().grid().comm());
1860 }
1861
1862
1863 template<typename TypeTag>
1864 typename BlackoilWellModel<TypeTag>::WellInterfacePtr
1865 BlackoilWellModel<TypeTag>::
1866 getWell(const std::string& well_name) const
1867 {
1868 // finding the iterator of the well in wells_ecl
1869 auto well = std::find_if(well_container_.begin(),
1870 well_container_.end(),
1871 [&well_name](const WellInterfacePtr& elem)->bool {
1872 return elem->name() == well_name;
1873 });
1874
1875 assert(well != well_container_.end());
1876
1877 return *well;
1878 }
1879
1880 template<typename TypeTag>
1881 bool
1882 BlackoilWellModel<TypeTag>::
1883 hasWell(const std::string& well_name) const
1884 {
1885 return std::any_of(well_container_.begin(), well_container_.end(),
1886 [&well_name](const WellInterfacePtr& elem) -> bool
1887 {
1888 return elem->name() == well_name;
1889 });
1890 }
1891
1892
1893
1894
1895 template <typename TypeTag>
1896 int
1897 BlackoilWellModel<TypeTag>::
1898 reportStepIndex() const
1899 {
1900 return std::max(this->ebosSimulator_.episodeIndex(), 0);
1901 }
1902
1903
1904
1905
1906
1907 template<typename TypeTag>
1908 void
1909 BlackoilWellModel<TypeTag>::
1910 calcRates(const int fipnum,
1911 const int pvtreg,
1912 std::vector<double>& resv_coeff)
1913 {
1914 rateConverter_->calcCoeff(fipnum, pvtreg, resv_coeff);
1915 }
1916
1917 template<typename TypeTag>
1918 void
1919 BlackoilWellModel<TypeTag>::
1920 calcInjRates(const int fipnum,
1921 const int pvtreg,
1922 std::vector<double>& resv_coeff)
1923 {
1924 rateConverter_->calcInjCoeff(fipnum, pvtreg, resv_coeff);
1925 }
1926
1927
1928 template <typename TypeTag>
1929 void
1930 BlackoilWellModel<TypeTag>::
1931 computeWellTemperature()
1932 {
1933 if (!has_energy_)
1934 return;
1935
1936 int np = numPhases();
1937 double cellInternalEnergy;
1938 double cellBinv;
1939 double cellDensity;
1940 double perfPhaseRate;
1941 const int nw = numLocalWells();
1942 for (auto wellID = 0*nw; wellID < nw; ++wellID) {
1943 const Well& well = wells_ecl_[wellID];
1944 if (well.isInjector())
1945 continue;
1946
1947 int connpos = 0;
1948 for (int i = 0; i < wellID; ++i) {
1949 connpos += well_perf_data_[i].size();
1950 }
1951 connpos *= np;
1952 std::array<double,2> weighted{0.0,0.0};
1953 auto& [weighted_temperature, total_weight] = weighted;
1954
1955 auto& well_info = local_parallel_well_info_[wellID].get();
1956 const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size());
1957 auto& ws = this->wellState().well(wellID);
1958 auto& perf_data = ws.perf_data;
1959 auto& perf_phase_rate = perf_data.phase_rates;
1960
1961 for (int perf = 0; perf < num_perf_this_well; ++perf) {
1962 const int cell_idx = well_perf_data_[wellID][perf].cell_index;
1963 const auto& intQuants = *(ebosSimulator_.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1964 const auto& fs = intQuants.fluidState();
1965
1966 // we on only have one temperature pr cell any phaseIdx will do
1967 double cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
1968
1969 double weight_factor = 0.0;
1970 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1971 {
1972 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1973 continue;
1974 }
1975 cellInternalEnergy = fs.enthalpy(phaseIdx).value() - fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
1976 cellBinv = fs.invB(phaseIdx).value();
1977 cellDensity = fs.density(phaseIdx).value();
1978 perfPhaseRate = perf_phase_rate[ perf*np + phaseIdx ];
1979 weight_factor += cellDensity * perfPhaseRate/cellBinv * cellInternalEnergy/cellTemperatures;
1980 }
1981 total_weight += weight_factor;
1982 weighted_temperature += weight_factor * cellTemperatures;
1983 }
1984 well_info.communication().sum(weighted.data(), 2);
1985 this->wellState().well(wellID).temperature = weighted_temperature/total_weight;
1986 }
1987 }
1988
1989
1990
1991 template <typename TypeTag>
1992 void
1993 BlackoilWellModel<TypeTag>::
1994 assignWellTracerRates(data::Wells& wsrpt) const
1995 {
1996 const auto & wellTracerRates = ebosSimulator_.problem().tracerModel().getWellTracerRates();
1997
1998 if (wellTracerRates.empty())
1999 return; // no tracers
2000
2001 for (const auto& wTR : wellTracerRates) {
2002 std::string wellName = wTR.first.first;
2003 auto xwPos = wsrpt.find(wellName);
2004 if (xwPos == wsrpt.end()) { // No well results.
2005 continue;
2006 }
2007 std::string tracerName = wTR.first.second;
2008 double rate = wTR.second;
2009 xwPos->second.rates.set(data::Rates::opt::tracer, rate, tracerName);
2010 }
2011 }
2012
2013
2014
2015
2016
2017} // namespace Opm
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:94
void calculateExplicitQuantities(DeferredLogger &deferred_logger) const
Calculating the explict quantities used in the well calculation.
Definition: BlackoilWellModel_impl.hpp:1445
Definition: DeferredLogger.hpp:57
void logMessages()
Log all messages to the OpmLog backends, and clear the message container.
Definition: DeferredLogger.cpp:85
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellState.hpp:56
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
Opm::DeferredLogger gatherDeferredLogger(const Opm::DeferredLogger &local_deferredlogger, Opm::Parallel::Communication)
Create a global log combining local logs.
Definition: gatherDeferredLogger.cpp:168
ConvergenceReport gatherConvergenceReport(const ConvergenceReport &local_report, Parallel::Communication mpi_communicator)
Create a global convergence report combining local (per-process) reports.
Definition: gatherConvergenceReport.cpp:205
PhaseUsage phaseUsageFromDeck(const EclipseState &eclipseState)
Looks at presence of WATER, OIL and GAS keywords in state object to determine active phases.
Definition: phaseUsageFromDeck.cpp:145