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