My Project
MultisegmentWell_impl.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #include <opm/simulators/wells/MSWellHelpers.hpp>
23 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>
25 #include <opm/common/OpmLog/OpmLog.hpp>
26 
27 #include <string>
28 #include <algorithm>
29 
30 #if HAVE_CUDA || HAVE_OPENCL
31 #include <opm/simulators/linalg/bda/WellContributions.hpp>
32 #endif
33 
34 namespace Opm
35 {
36 
37 
38  template <typename TypeTag>
39  MultisegmentWell<TypeTag>::
40  MultisegmentWell(const Well& well,
41  const ParallelWellInfo& pw_info,
42  const int time_step,
43  const ModelParameters& param,
44  const RateConverterType& rate_converter,
45  const int pvtRegionIdx,
46  const int num_components,
47  const int num_phases,
48  const int index_of_well,
49  const std::vector<PerforationData>& perf_data)
50  : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
51  , MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
52  , segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
53  {
54  // not handling solvent or polymer for now with multisegment well
55  if constexpr (has_solvent) {
56  OPM_THROW(std::runtime_error, "solvent is not supported by multisegment well yet");
57  }
58 
59  if constexpr (has_polymer) {
60  OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet");
61  }
62 
63  if constexpr (Base::has_energy) {
64  OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet");
65  }
66 
67  if constexpr (Base::has_foam) {
68  OPM_THROW(std::runtime_error, "foam is not supported by multisegment well yet");
69  }
70 
71  if constexpr (Base::has_brine) {
72  OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet");
73  }
74  }
75 
76 
77 
78 
79 
80  template <typename TypeTag>
81  void
82  MultisegmentWell<TypeTag>::
83  init(const PhaseUsage* phase_usage_arg,
84  const std::vector<double>& depth_arg,
85  const double gravity_arg,
86  const int num_cells,
87  const std::vector< Scalar >& B_avg)
88  {
89  Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg);
90 
91  // TODO: for StandardWell, we need to update the perf depth here using depth_arg.
92  // for MultisegmentWell, it is much more complicated.
93  // It can be specified directly, it can be calculated from the segment depth,
94  // it can also use the cell center, which is the same for StandardWell.
95  // For the last case, should we update the depth with the depth_arg? For the
96  // future, it can be a source of wrong result with Multisegment well.
97  // An indicator from the opm-parser should indicate what kind of depth we should use here.
98 
99  // \Note: we do not update the depth here. And it looks like for now, we only have the option to use
100  // specified perforation depth
101  this->initMatrixAndVectors(num_cells);
102 
103  // calcuate the depth difference between the perforations and the perforated grid block
104  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
105  const int cell_idx = this->well_cells_[perf];
106  this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf];
107  }
108  }
109 
110 
111 
112 
113 
114  template <typename TypeTag>
115  void
116  MultisegmentWell<TypeTag>::
117  initPrimaryVariablesEvaluation() const
118  {
119  this->MSWEval::initPrimaryVariablesEvaluation();
120  }
121 
122 
123 
124 
125 
126  template <typename TypeTag>
127  void
128  MultisegmentWell<TypeTag>::
129  updatePrimaryVariables(const WellState& well_state, DeferredLogger& /* deferred_logger */) const
130  {
131  this->MSWEval::updatePrimaryVariables(well_state);
132  }
133 
134 
135 
136 
137 
138 
139  template <typename TypeTag>
140  void
142  updateWellStateWithTarget(const Simulator& ebos_simulator,
143  const GroupState& group_state,
144  WellState& well_state,
145  DeferredLogger& deferred_logger) const
146  {
147  Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
148  // scale segment rates based on the wellRates
149  // and segment pressure based on bhp
150  this->scaleSegmentRatesWithWellRates(well_state);
151  this->scaleSegmentPressuresWithBhp(well_state);
152  }
153 
154 
155 
156 
157 
158  template <typename TypeTag>
161  getWellConvergence(const WellState& well_state,
162  const std::vector<double>& B_avg,
163  DeferredLogger& deferred_logger,
164  const bool relax_tolerance) const
165  {
166  return this->MSWEval::getWellConvergence(well_state,
167  B_avg,
168  deferred_logger,
169  this->param_.max_residual_allowed_,
170  this->param_.tolerance_wells_,
171  this->param_.relaxed_tolerance_flow_well_,
172  this->param_.tolerance_pressure_ms_wells_,
173  this->param_.relaxed_tolerance_pressure_ms_well_,
174  relax_tolerance);
175  }
176 
177 
178 
179 
180 
181  template <typename TypeTag>
182  void
184  apply(const BVector& x, BVector& Ax) const
185  {
186  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
187 
188  if ( this->param_.matrix_add_well_contributions_ )
189  {
190  // Contributions are already in the matrix itself
191  return;
192  }
193  BVectorWell Bx(this->duneB_.N());
194 
195  this->duneB_.mv(x, Bx);
196 
197  // invDBx = duneD^-1 * Bx_
198  const BVectorWell invDBx = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, Bx);
199 
200  // Ax = Ax - duneC_^T * invDBx
201  this->duneC_.mmtv(invDBx,Ax);
202  }
203 
204 
205 
206 
207 
208  template <typename TypeTag>
209  void
211  apply(BVector& r) const
212  {
213  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
214 
215  // invDrw_ = duneD^-1 * resWell_
216  const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
217  // r = r - duneC_^T * invDrw
218  this->duneC_.mmtv(invDrw, r);
219  }
220 
221 
222 
223  template <typename TypeTag>
224  void
227  WellState& well_state,
228  DeferredLogger& deferred_logger) const
229  {
230  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
231 
232  BVectorWell xw(1);
233  this->recoverSolutionWell(x, xw);
234  updateWellState(xw, well_state, deferred_logger);
235  }
236 
237 
238 
239 
240 
241  template <typename TypeTag>
242  void
244  computeWellPotentials(const Simulator& ebosSimulator,
245  const WellState& well_state,
246  std::vector<double>& well_potentials,
247  DeferredLogger& deferred_logger)
248  {
249  const int np = this->number_of_phases_;
250  well_potentials.resize(np, 0.0);
251 
252  // Stopped wells have zero potential.
253  if (this->wellIsStopped()) {
254  return;
255  }
256 
257  // If the well is pressure controlled the potential equals the rate.
258  bool thp_controlled_well = false;
259  bool bhp_controlled_well = false;
260  const auto& ws = well_state.well(this->index_of_well_);
261  if (this->isInjector()) {
262  const Well::InjectorCMode& current = ws.injection_cmode;
263  if (current == Well::InjectorCMode::THP) {
264  thp_controlled_well = true;
265  }
266  if (current == Well::InjectorCMode::BHP) {
267  bhp_controlled_well = true;
268  }
269  } else {
270  const Well::ProducerCMode& current = ws.production_cmode;
271  if (current == Well::ProducerCMode::THP) {
272  thp_controlled_well = true;
273  }
274  if (current == Well::ProducerCMode::BHP) {
275  bhp_controlled_well = true;
276  }
277  }
278  if (thp_controlled_well || bhp_controlled_well) {
279 
280  double total_rate = 0.0;
281  for (int phase = 0; phase < np; ++phase){
282  total_rate += ws.surface_rates[phase];
283  }
284  // for pressure controlled wells the well rates are the potentials
285  // if the rates are trivial we are most probably looking at the newly
286  // opened well and we therefore make the affort of computing the potentials anyway.
287  if (std::abs(total_rate) > 0) {
288  for (int phase = 0; phase < np; ++phase){
289  well_potentials[phase] = ws.surface_rates[phase];
290  }
291  return;
292  }
293  }
294 
295  debug_cost_counter_ = 0;
296  // does the well have a THP related constraint?
297  const auto& summaryState = ebosSimulator.vanguard().summaryState();
298  if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
299  computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
300  } else {
301  well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
302  }
303  deferred_logger.debug("Cost in iterations of finding well potential for well "
304  + this->name() + ": " + std::to_string(debug_cost_counter_));
305  }
306 
307 
308 
309 
310  template<typename TypeTag>
311  void
313  computeWellRatesAtBhpLimit(const Simulator& ebosSimulator,
314  std::vector<double>& well_flux,
315  DeferredLogger& deferred_logger) const
316  {
317  if (this->well_ecl_.isInjector()) {
318  const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState());
319  computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
320  } else {
321  const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState());
322  computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger);
323  }
324  }
325 
326  template<typename TypeTag>
327  void
328  MultisegmentWell<TypeTag>::
329  computeWellRatesWithBhp(const Simulator& ebosSimulator,
330  const Scalar& bhp,
331  std::vector<double>& well_flux,
332  DeferredLogger& deferred_logger) const
333  {
334 
335  const int np = this->number_of_phases_;
336 
337  well_flux.resize(np, 0.0);
338  const bool allow_cf = this->getAllowCrossFlow();
339  const int nseg = this->numberOfSegments();
340  const WellState& well_state = ebosSimulator.problem().wellModel().wellState();
341  const auto& ws = well_state.well(this->indexOfWell());
342  auto segments_copy = ws.segments;
343  segments_copy.scale_pressure(bhp);
344  const auto& segment_pressure = segments_copy.pressure;
345  for (int seg = 0; seg < nseg; ++seg) {
346  for (const int perf : this->segment_perforations_[seg]) {
347  const int cell_idx = this->well_cells_[perf];
348  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
349  // flux for each perforation
350  std::vector<Scalar> mob(this->num_components_, 0.);
351  getMobilityScalar(ebosSimulator, perf, mob);
352  double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
353  const double Tw = this->well_index_[perf] * trans_mult;
354 
355  const Scalar seg_pressure = segment_pressure[seg];
356  std::vector<Scalar> cq_s(this->num_components_, 0.);
357  computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure,
358  allow_cf, cq_s, deferred_logger);
359 
360  for(int p = 0; p < np; ++p) {
361  well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
362  }
363  }
364  }
365  this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
366  }
367 
368 
369  template<typename TypeTag>
370  void
371  MultisegmentWell<TypeTag>::
372  computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
373  const Scalar& bhp,
374  std::vector<double>& well_flux,
375  DeferredLogger& deferred_logger) const
376  {
377  // creating a copy of the well itself, to avoid messing up the explicit informations
378  // during this copy, the only information not copied properly is the well controls
379  MultisegmentWell<TypeTag> well_copy(*this);
380  well_copy.debug_cost_counter_ = 0;
381 
382  // store a copy of the well state, we don't want to update the real well state
383  WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
384  const auto& group_state = ebosSimulator.problem().wellModel().groupState();
385  auto& ws = well_state_copy.well(this->index_of_well_);
386 
387  // Get the current controls.
388  const auto& summary_state = ebosSimulator.vanguard().summaryState();
389  auto inj_controls = well_copy.well_ecl_.isInjector()
390  ? well_copy.well_ecl_.injectionControls(summary_state)
391  : Well::InjectionControls(0);
392  auto prod_controls = well_copy.well_ecl_.isProducer()
393  ? well_copy.well_ecl_.productionControls(summary_state) :
394  Well::ProductionControls(0);
395 
396  // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
397  if (well_copy.well_ecl_.isInjector()) {
398  inj_controls.bhp_limit = bhp;
399  ws.injection_cmode = Well::InjectorCMode::BHP;
400  } else {
401  prod_controls.bhp_limit = bhp;
402  ws.production_cmode = Well::ProducerCMode::BHP;
403  }
404  ws.bhp = bhp;
405  well_copy.scaleSegmentPressuresWithBhp(well_state_copy);
406 
407  // initialized the well rates with the potentials i.e. the well rates based on bhp
408  const int np = this->number_of_phases_;
409  const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
410  for (int phase = 0; phase < np; ++phase){
411  ws.surface_rates[phase] = sign * ws.well_potentials[phase];
412  }
413  well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
414 
415  well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
416  const double dt = ebosSimulator.timeStepSize();
417  // iterate to get a solution at the given bhp.
418  well_copy.iterateWellEqWithControl(ebosSimulator, dt, inj_controls, prod_controls, well_state_copy, group_state,
419  deferred_logger);
420 
421  // compute the potential and store in the flux vector.
422  well_flux.clear();
423  well_flux.resize(np, 0.0);
424  for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
425  const EvalWell rate = well_copy.getQs(compIdx);
426  well_flux[this->ebosCompIdxToFlowCompIdx(compIdx)] = rate.value();
427  }
428  debug_cost_counter_ += well_copy.debug_cost_counter_;
429  }
430 
431 
432 
433  template<typename TypeTag>
434  std::vector<double>
435  MultisegmentWell<TypeTag>::
436  computeWellPotentialWithTHP(const Simulator& ebos_simulator,
437  DeferredLogger& deferred_logger) const
438  {
439  std::vector<double> potentials(this->number_of_phases_, 0.0);
440  const auto& summary_state = ebos_simulator.vanguard().summaryState();
441 
442  const auto& well = this->well_ecl_;
443  if (well.isInjector()){
444  auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
445  if (bhp_at_thp_limit) {
446  const auto& controls = well.injectionControls(summary_state);
447  const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
448  computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
449  deferred_logger.debug("Converged thp based potential calculation for well "
450  + this->name() + ", at bhp = " + std::to_string(bhp));
451  } else {
452  deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
453  "Failed in getting converged thp based potential calculation for well "
454  + this->name() + ". Instead the bhp based value is used");
455  const auto& controls = well.injectionControls(summary_state);
456  const double bhp = controls.bhp_limit;
457  computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
458  }
459  } else {
460  auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
461  if (bhp_at_thp_limit) {
462  const auto& controls = well.productionControls(summary_state);
463  const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
464  computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
465  deferred_logger.debug("Converged thp based potential calculation for well "
466  + this->name() + ", at bhp = " + std::to_string(bhp));
467  } else {
468  deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
469  "Failed in getting converged thp based potential calculation for well "
470  + this->name() + ". Instead the bhp based value is used");
471  const auto& controls = well.productionControls(summary_state);
472  const double bhp = controls.bhp_limit;
473  computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
474  }
475  }
476 
477  return potentials;
478  }
479 
480 
481 
482  template <typename TypeTag>
483  void
484  MultisegmentWell<TypeTag>::
485  solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
486  {
487  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
488 
489  // We assemble the well equations, then we check the convergence,
490  // which is why we do not put the assembleWellEq here.
491  const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
492 
493  updateWellState(dx_well, well_state, deferred_logger);
494  }
495 
496 
497 
498 
499 
500  template <typename TypeTag>
501  void
502  MultisegmentWell<TypeTag>::
503  computePerfCellPressDiffs(const Simulator& ebosSimulator)
504  {
505  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
506 
507  std::vector<double> kr(this->number_of_phases_, 0.0);
508  std::vector<double> density(this->number_of_phases_, 0.0);
509 
510  const int cell_idx = this->well_cells_[perf];
511  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
512  const auto& fs = intQuants.fluidState();
513 
514  double sum_kr = 0.;
515 
516  const PhaseUsage& pu = this->phaseUsage();
517  if (pu.phase_used[Water]) {
518  const int water_pos = pu.phase_pos[Water];
519  kr[water_pos] = intQuants.relativePermeability(FluidSystem::waterPhaseIdx).value();
520  sum_kr += kr[water_pos];
521  density[water_pos] = fs.density(FluidSystem::waterPhaseIdx).value();
522  }
523 
524  if (pu.phase_used[Oil]) {
525  const int oil_pos = pu.phase_pos[Oil];
526  kr[oil_pos] = intQuants.relativePermeability(FluidSystem::oilPhaseIdx).value();
527  sum_kr += kr[oil_pos];
528  density[oil_pos] = fs.density(FluidSystem::oilPhaseIdx).value();
529  }
530 
531  if (pu.phase_used[Gas]) {
532  const int gas_pos = pu.phase_pos[Gas];
533  kr[gas_pos] = intQuants.relativePermeability(FluidSystem::gasPhaseIdx).value();
534  sum_kr += kr[gas_pos];
535  density[gas_pos] = fs.density(FluidSystem::gasPhaseIdx).value();
536  }
537 
538  assert(sum_kr != 0.);
539 
540  // calculate the average density
541  double average_density = 0.;
542  for (int p = 0; p < this->number_of_phases_; ++p) {
543  average_density += kr[p] * density[p];
544  }
545  average_density /= sum_kr;
546 
547  this->cell_perforation_pressure_diffs_[perf] = this->gravity_ * average_density * this->cell_perforation_depth_diffs_[perf];
548  }
549  }
550 
551 
552 
553 
554 
555  template <typename TypeTag>
556  void
557  MultisegmentWell<TypeTag>::
558  computeInitialSegmentFluids(const Simulator& ebos_simulator)
559  {
560  for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
561  // TODO: trying to reduce the times for the surfaceVolumeFraction calculation
562  const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value();
563  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
564  segment_fluid_initial_[seg][comp_idx] = surface_volume * this->surfaceVolumeFraction(seg, comp_idx).value();
565  }
566  }
567  }
568 
569 
570 
571 
572 
573  template <typename TypeTag>
574  void
575  MultisegmentWell<TypeTag>::
576  updateWellState(const BVectorWell& dwells,
577  WellState& well_state,
578  DeferredLogger& deferred_logger,
579  const double relaxation_factor) const
580  {
581  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
582 
583  const double dFLimit = this->param_.dwell_fraction_max_;
584  const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
585  this->MSWEval::updateWellState(dwells,
586  relaxation_factor,
587  dFLimit,
588  max_pressure_change);
589 
590  this->updateWellStateFromPrimaryVariables(well_state, getRefDensity(), deferred_logger);
591  Base::calculateReservoirRates(well_state.well(this->index_of_well_));
592  }
593 
594 
595 
596 
597 
598  template <typename TypeTag>
599  void
600  MultisegmentWell<TypeTag>::
601  calculateExplicitQuantities(const Simulator& ebosSimulator,
602  const WellState& well_state,
603  DeferredLogger& deferred_logger)
604  {
605  updatePrimaryVariables(well_state, deferred_logger);
606  initPrimaryVariablesEvaluation();
607  computePerfCellPressDiffs(ebosSimulator);
608  computeInitialSegmentFluids(ebosSimulator);
609  }
610 
611 
612 
613 
614 
615  template<typename TypeTag>
616  void
617  MultisegmentWell<TypeTag>::
618  updateProductivityIndex(const Simulator& ebosSimulator,
619  const WellProdIndexCalculator& wellPICalc,
620  WellState& well_state,
621  DeferredLogger& deferred_logger) const
622  {
623  auto fluidState = [&ebosSimulator, this](const int perf)
624  {
625  const auto cell_idx = this->well_cells_[perf];
626  return ebosSimulator.model()
627  .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
628  };
629 
630  const int np = this->number_of_phases_;
631  auto setToZero = [np](double* x) -> void
632  {
633  std::fill_n(x, np, 0.0);
634  };
635 
636  auto addVector = [np](const double* src, double* dest) -> void
637  {
638  std::transform(src, src + np, dest, dest, std::plus<>{});
639  };
640 
641  auto& ws = well_state.well(this->index_of_well_);
642  auto& perf_data = ws.perf_data;
643  auto* connPI = perf_data.prod_index.data();
644  auto* wellPI = ws.productivity_index.data();
645 
646  setToZero(wellPI);
647 
648  const auto preferred_phase = this->well_ecl_.getPreferredPhase();
649  auto subsetPerfID = 0;
650 
651  for ( const auto& perf : *this->perf_data_){
652  auto allPerfID = perf.ecl_index;
653 
654  auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
655  {
656  return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
657  };
658 
659  std::vector<Scalar> mob(this->num_components_, 0.0);
660  getMobilityScalar(ebosSimulator, static_cast<int>(subsetPerfID), mob);
661 
662  const auto& fs = fluidState(subsetPerfID);
663  setToZero(connPI);
664 
665  if (this->isInjector()) {
666  this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
667  mob, connPI, deferred_logger);
668  }
669  else { // Production or zero flow rate
670  this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
671  }
672 
673  addVector(connPI, wellPI);
674 
675  ++subsetPerfID;
676  connPI += np;
677  }
678 
679  assert (static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
680  "Internal logic error in processing connections for PI/II");
681  }
682 
683 
684 
685 
686 
687  template<typename TypeTag>
688  void
689  MultisegmentWell<TypeTag>::
690  addWellContributions(SparseMatrixAdapter& jacobian) const
691  {
692  const auto invDuneD = mswellhelpers::invertWithUMFPack<DiagMatWell, BVectorWell>(this->duneD_, this->duneDSolver_);
693 
694  // We need to change matrix A as follows
695  // A -= C^T D^-1 B
696  // D is a (nseg x nseg) block matrix with (4 x 4) blocks.
697  // B and C are (nseg x ncells) block matrices with (4 x 4 blocks).
698  // They have nonzeros at (i, j) only if this well has a
699  // perforation at cell j connected to segment i. The code
700  // assumes that no cell is connected to more than one segment,
701  // i.e. the columns of B/C have no more than one nonzero.
702  for (size_t rowC = 0; rowC < this->duneC_.N(); ++rowC) {
703  for (auto colC = this->duneC_[rowC].begin(), endC = this->duneC_[rowC].end(); colC != endC; ++colC) {
704  const auto row_index = colC.index();
705  for (size_t rowB = 0; rowB < this->duneB_.N(); ++rowB) {
706  for (auto colB = this->duneB_[rowB].begin(), endB = this->duneB_[rowB].end(); colB != endB; ++colB) {
707  const auto col_index = colB.index();
708  OffDiagMatrixBlockWellType tmp1;
709  Detail::multMatrixImpl(invDuneD[rowC][rowB], (*colB), tmp1, std::true_type());
710  typename SparseMatrixAdapter::MatrixBlock tmp2;
711  Detail::multMatrixTransposedImpl((*colC), tmp1, tmp2, std::false_type());
712  jacobian.addToBlock(row_index, col_index, tmp2);
713  }
714  }
715  }
716  }
717  }
718 
719 
720 
721  template<typename TypeTag>
722  template<class Value>
723  void
724  MultisegmentWell<TypeTag>::
725  computePerfRate(const Value& pressure_cell,
726  const Value& rs,
727  const Value& rv,
728  const std::vector<Value>& b_perfcells,
729  const std::vector<Value>& mob_perfcells,
730  const double Tw,
731  const int perf,
732  const Value& segment_pressure,
733  const Value& segment_density,
734  const bool& allow_cf,
735  const std::vector<Value>& cmix_s,
736  std::vector<Value>& cq_s,
737  Value& perf_press,
738  double& perf_dis_gas_rate,
739  double& perf_vap_oil_rate,
740  DeferredLogger& deferred_logger) const
741  {
742  // pressure difference between the segment and the perforation
743  const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
744  // pressure difference between the perforation and the grid cell
745  const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
746 
747  perf_press = pressure_cell - cell_perf_press_diff;
748  // Pressure drawdown (also used to determine direction of flow)
749  // TODO: not 100% sure about the sign of the seg_perf_press_diff
750  const Value drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
751 
752  // producing perforations
753  if ( drawdown > 0.0) {
754  // Do nothing is crossflow is not allowed
755  if (!allow_cf && this->isInjector()) {
756  return;
757  }
758 
759  // compute component volumetric rates at standard conditions
760  for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
761  const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
762  cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
763  }
764 
765  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
766  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
767  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
768  const Value cq_s_oil = cq_s[oilCompIdx];
769  const Value cq_s_gas = cq_s[gasCompIdx];
770  cq_s[gasCompIdx] += rs * cq_s_oil;
771  cq_s[oilCompIdx] += rv * cq_s_gas;
772  }
773  } else { // injecting perforations
774  // Do nothing if crossflow is not allowed
775  if (!allow_cf && this->isProducer()) {
776  return;
777  }
778 
779  // for injecting perforations, we use total mobility
780  Value total_mob = mob_perfcells[0];
781  for (int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
782  total_mob += mob_perfcells[comp_idx];
783  }
784 
785  // injection perforations total volume rates
786  const Value cqt_i = - Tw * (total_mob * drawdown);
787 
788  // compute volume ratio between connection and at standard conditions
789  Value volume_ratio = 0.0;
790  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
791  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
792  volume_ratio += cmix_s[waterCompIdx] / b_perfcells[waterCompIdx];
793  }
794 
795  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
796  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
797  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
798 
799  // Incorporate RS/RV factors if both oil and gas active
800  // TODO: not sure we use rs rv from the perforation cells when handling injecting perforations
801  // basically, for injecting perforations, the wellbore is the upstreaming side.
802  const Value d = 1.0 - rv * rs;
803 
804  if (getValue(d) == 0.0) {
805  OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << this->name()
806  << " during flux calculation"
807  << " with rs " << rs << " and rv " << rv, deferred_logger);
808  }
809 
810  const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
811  volume_ratio += tmp_oil / b_perfcells[oilCompIdx];
812 
813  const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
814  volume_ratio += tmp_gas / b_perfcells[gasCompIdx];
815  } else { // not having gas and oil at the same time
816  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
817  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
818  volume_ratio += cmix_s[oilCompIdx] / b_perfcells[oilCompIdx];
819  }
820  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
821  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
822  volume_ratio += cmix_s[gasCompIdx] / b_perfcells[gasCompIdx];
823  }
824  }
825  // injecting connections total volumerates at standard conditions
826  Value cqt_is = cqt_i / volume_ratio;
827  for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
828  cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
829  }
830  } // end for injection perforations
831 
832  // calculating the perforation solution gas rate and solution oil rates
833  if (this->isProducer()) {
834  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
835  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
836  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
837  // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
838  // s means standard condition, r means reservoir condition
839  // q_os = q_or * b_o + rv * q_gr * b_g
840  // q_gs = q_gr * g_g + rs * q_or * b_o
841  // d = 1.0 - rs * rv
842  // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
843  // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
844 
845  const double d = 1.0 - getValue(rv) * getValue(rs);
846  // vaporized oil into gas
847  // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
848  perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
849  // dissolved of gas in oil
850  // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
851  perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
852  }
853  }
854  }
855 
856  template <typename TypeTag>
857  void
858  MultisegmentWell<TypeTag>::
859  computePerfRateEval(const IntensiveQuantities& int_quants,
860  const std::vector<EvalWell>& mob_perfcells,
861  const double Tw,
862  const int seg,
863  const int perf,
864  const EvalWell& segment_pressure,
865  const bool& allow_cf,
866  std::vector<EvalWell>& cq_s,
867  EvalWell& perf_press,
868  double& perf_dis_gas_rate,
869  double& perf_vap_oil_rate,
870  DeferredLogger& deferred_logger) const
871 
872  {
873  const auto& fs = int_quants.fluidState();
874 
875  const EvalWell pressure_cell = this->extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
876  const EvalWell rs = this->extendEval(fs.Rs());
877  const EvalWell rv = this->extendEval(fs.Rv());
878 
879  // not using number_of_phases_ because of solvent
880  std::vector<EvalWell> b_perfcells(this->num_components_, 0.0);
881 
882  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
883  if (!FluidSystem::phaseIsActive(phaseIdx)) {
884  continue;
885  }
886 
887  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
888  b_perfcells[compIdx] = this->extendEval(fs.invB(phaseIdx));
889  }
890 
891  std::vector<EvalWell> cmix_s(this->numComponents(), 0.0);
892  for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
893  cmix_s[comp_idx] = this->surfaceVolumeFraction(seg, comp_idx);
894  }
895 
896  this->computePerfRate(pressure_cell,
897  rs,
898  rv,
899  b_perfcells,
900  mob_perfcells,
901  Tw,
902  perf,
903  segment_pressure,
904  this->segment_densities_[seg],
905  allow_cf,
906  cmix_s,
907  cq_s,
908  perf_press,
909  perf_dis_gas_rate,
910  perf_vap_oil_rate,
911  deferred_logger);
912  }
913 
914 
915 
916  template <typename TypeTag>
917  void
918  MultisegmentWell<TypeTag>::
919  computePerfRateScalar(const IntensiveQuantities& int_quants,
920  const std::vector<Scalar>& mob_perfcells,
921  const double Tw,
922  const int seg,
923  const int perf,
924  const Scalar& segment_pressure,
925  const bool& allow_cf,
926  std::vector<Scalar>& cq_s,
927  DeferredLogger& deferred_logger) const
928 
929  {
930  const auto& fs = int_quants.fluidState();
931 
932  const Scalar pressure_cell = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
933  const Scalar rs = getValue(fs.Rs());
934  const Scalar rv = getValue(fs.Rv());
935 
936  // not using number_of_phases_ because of solvent
937  std::vector<Scalar> b_perfcells(this->num_components_, 0.0);
938 
939  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
940  if (!FluidSystem::phaseIsActive(phaseIdx)) {
941  continue;
942  }
943 
944  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
945  b_perfcells[compIdx] = getValue(fs.invB(phaseIdx));
946  }
947 
948  std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
949  for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
950  cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx));
951  }
952 
953  Scalar perf_dis_gas_rate = 0.0;
954  Scalar perf_vap_oil_rate = 0.0;
955  Scalar perf_press = 0.0;
956 
957  this->computePerfRate(pressure_cell,
958  rs,
959  rv,
960  b_perfcells,
961  mob_perfcells,
962  Tw,
963  perf,
964  segment_pressure,
965  getValue(this->segment_densities_[seg]),
966  allow_cf,
967  cmix_s,
968  cq_s,
969  perf_press,
970  perf_dis_gas_rate,
971  perf_vap_oil_rate,
972  deferred_logger);
973  }
974 
975  template <typename TypeTag>
976  void
977  MultisegmentWell<TypeTag>::
978  computeSegmentFluidProperties(const Simulator& ebosSimulator)
979  {
980  // TODO: the concept of phases and components are rather confusing in this function.
981  // needs to be addressed sooner or later.
982 
983  // get the temperature for later use. It is only useful when we are not handling
984  // thermal related simulation
985  // basically, it is a single value for all the segments
986 
987  EvalWell temperature;
988  EvalWell saltConcentration;
989  // not sure how to handle the pvt region related to segment
990  // for the current approach, we use the pvt region of the first perforated cell
991  // although there are some text indicating using the pvt region of the lowest
992  // perforated cell
993  // TODO: later to investigate how to handle the pvt region
994  int pvt_region_index;
995  {
996  // using the first perforated cell
997  const int cell_idx = this->well_cells_[0];
998  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
999  const auto& fs = intQuants.fluidState();
1000  temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1001  saltConcentration = this->extendEval(fs.saltConcentration());
1002  pvt_region_index = fs.pvtRegionIndex();
1003  }
1004 
1005  this->MSWEval::computeSegmentFluidProperties(temperature,
1006  saltConcentration,
1007  pvt_region_index);
1008  }
1009 
1010 
1011 
1012 
1013 
1014  template <typename TypeTag>
1015  void
1016  MultisegmentWell<TypeTag>::
1017  getMobilityEval(const Simulator& ebosSimulator,
1018  const int perf,
1019  std::vector<EvalWell>& mob) const
1020  {
1021  // TODO: most of this function, if not the whole function, can be moved to the base class
1022  const int cell_idx = this->well_cells_[perf];
1023  assert (int(mob.size()) == this->num_components_);
1024  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1025  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1026 
1027  // either use mobility of the perforation cell or calcualte its own
1028  // based on passing the saturation table index
1029  const int satid = this->saturation_table_number_[perf] - 1;
1030  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1031  if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
1032 
1033  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1034  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1035  continue;
1036  }
1037 
1038  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1039  mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
1040  }
1041  // if (has_solvent) {
1042  // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1043  // }
1044  } else {
1045 
1046  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1047  Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
1048  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1049 
1050  // reset the satnumvalue back to original
1051  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1052 
1053  // compute the mobility
1054  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1055  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1056  continue;
1057  }
1058 
1059  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1060  mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
1061  }
1062  }
1063  }
1064 
1065 
1066  template <typename TypeTag>
1067  void
1068  MultisegmentWell<TypeTag>::
1069  getMobilityScalar(const Simulator& ebosSimulator,
1070  const int perf,
1071  std::vector<Scalar>& mob) const
1072  {
1073  // TODO: most of this function, if not the whole function, can be moved to the base class
1074  const int cell_idx = this->well_cells_[perf];
1075  assert (int(mob.size()) == this->num_components_);
1076  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1077  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
1078 
1079  // either use mobility of the perforation cell or calcualte its own
1080  // based on passing the saturation table index
1081  const int satid = this->saturation_table_number_[perf] - 1;
1082  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
1083  if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
1084 
1085  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1086  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1087  continue;
1088  }
1089 
1090  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1091  mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
1092  }
1093  // if (has_solvent) {
1094  // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
1095  // }
1096  } else {
1097 
1098  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
1099  Scalar relativePerms[3] = { 0.0, 0.0, 0.0 };
1100  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
1101 
1102  // reset the satnumvalue back to original
1103  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
1104 
1105  // compute the mobility
1106  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1107  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1108  continue;
1109  }
1110 
1111  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1112  mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
1113  }
1114  }
1115  }
1116 
1117 
1118 
1119 
1120  template<typename TypeTag>
1121  double
1122  MultisegmentWell<TypeTag>::
1123  getRefDensity() const
1124  {
1125  return this->segment_densities_[0].value();
1126  }
1127 
1128  template<typename TypeTag>
1129  void
1130  MultisegmentWell<TypeTag>::
1131  checkOperabilityUnderBHPLimitProducer(const WellState& /*well_state*/, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1132  {
1133  const auto& summaryState = ebos_simulator.vanguard().summaryState();
1134  const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1135  // Crude but works: default is one atmosphere.
1136  // TODO: a better way to detect whether the BHP is defaulted or not
1137  const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1138  if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1139  // if the BHP limit is not defaulted or the well does not have a THP limit
1140  // we need to check the BHP limit
1141 
1142  double temp = 0;
1143  for (int p = 0; p < this->number_of_phases_; ++p) {
1144  temp += this->ipr_a_[p] - this->ipr_b_[p] * bhp_limit;
1145  }
1146  if (temp < 0.) {
1147  this->operability_status_.operable_under_only_bhp_limit = false;
1148  }
1149 
1150  // checking whether running under BHP limit will violate THP limit
1151  if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1152  // option 1: calculate well rates based on the BHP limit.
1153  // option 2: stick with the above IPR curve
1154  // we use IPR here
1155  std::vector<double> well_rates_bhp_limit;
1156  computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1157 
1158  const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger);
1159 
1160  const double thp_limit = this->getTHPConstraint(summaryState);
1161 
1162  if (thp < thp_limit) {
1163  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1164  }
1165  }
1166  } else {
1167  // defaulted BHP and there is a THP constraint
1168  // default BHP limit is about 1 atm.
1169  // when applied the hydrostatic pressure correction dp,
1170  // most likely we get a negative value (bhp + dp)to search in the VFP table,
1171  // which is not desirable.
1172  // we assume we can operate under defaulted BHP limit and will violate the THP limit
1173  // when operating under defaulted BHP limit.
1174  this->operability_status_.operable_under_only_bhp_limit = true;
1175  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1176  }
1177  }
1178 
1179 
1180 
1181  template<typename TypeTag>
1182  void
1183  MultisegmentWell<TypeTag>::
1184  updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const
1185  {
1186  // TODO: not handling solvent related here for now
1187 
1188  // TODO: it only handles the producers for now
1189  // the formular for the injectors are not formulated yet
1190  if (this->isInjector()) {
1191  return;
1192  }
1193 
1194  // initialize all the values to be zero to begin with
1195  std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1196  std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1197 
1198  const int nseg = this->numberOfSegments();
1199  double seg_bhp_press_diff = 0;
1200  double ref_depth = this->ref_depth_;
1201  for (int seg = 0; seg < nseg; ++seg) {
1202  // calculating the perforation rate for each perforation that belongs to this segment
1203  const double segment_depth = this->segmentSet()[seg].depth();
1204  const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, this->segment_densities_[seg].value(), this->gravity_);
1205  ref_depth = segment_depth;
1206  seg_bhp_press_diff += dp;
1207  for (const int perf : this->segment_perforations_[seg]) {
1208  std::vector<Scalar> mob(this->num_components_, 0.0);
1209 
1210  // TODO: mabye we should store the mobility somewhere, so that we only need to calculate it one per iteration
1211  getMobilityScalar(ebos_simulator, perf, mob);
1212 
1213  const int cell_idx = this->well_cells_[perf];
1214  const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1215  const auto& fs = int_quantities.fluidState();
1216  // the pressure of the reservoir grid block the well connection is in
1217  // pressure difference between the segment and the perforation
1218  const double perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
1219  // pressure difference between the perforation and the grid cell
1220  const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1221  const double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
1222 
1223  // calculating the b for the connection
1224  std::vector<double> b_perf(this->num_components_);
1225  for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1226  if (!FluidSystem::phaseIsActive(phase)) {
1227  continue;
1228  }
1229  const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1230  b_perf[comp_idx] = fs.invB(phase).value();
1231  }
1232 
1233  // the pressure difference between the connection and BHP
1234  const double h_perf = cell_perf_press_diff + perf_seg_press_diff + seg_bhp_press_diff;
1235  const double pressure_diff = pressure_cell - h_perf;
1236 
1237  // Let us add a check, since the pressure is calculated based on zero value BHP
1238  // it should not be negative anyway. If it is negative, we might need to re-formulate
1239  // to taking into consideration the crossflow here.
1240  if (pressure_diff <= 0.) {
1241  deferred_logger.warning("NON_POSITIVE_DRAWDOWN_IPR",
1242  "non-positive drawdown found when updateIPR for well " + this->name());
1243  }
1244 
1245  // the well index associated with the connection
1246  const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1247 
1248  // TODO: there might be some indices related problems here
1249  // phases vs components
1250  // ipr values for the perforation
1251  std::vector<double> ipr_a_perf(this->ipr_a_.size());
1252  std::vector<double> ipr_b_perf(this->ipr_b_.size());
1253  for (int p = 0; p < this->number_of_phases_; ++p) {
1254  const double tw_mob = tw_perf * mob[p] * b_perf[p];
1255  ipr_a_perf[p] += tw_mob * pressure_diff;
1256  ipr_b_perf[p] += tw_mob;
1257  }
1258 
1259  // we need to handle the rs and rv when both oil and gas are present
1260  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1261  const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1262  const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1263  const double rs = (fs.Rs()).value();
1264  const double rv = (fs.Rv()).value();
1265 
1266  const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1267  const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1268 
1269  ipr_a_perf[gas_comp_idx] += dis_gas_a;
1270  ipr_a_perf[oil_comp_idx] += vap_oil_a;
1271 
1272  const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1273  const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1274 
1275  ipr_b_perf[gas_comp_idx] += dis_gas_b;
1276  ipr_b_perf[oil_comp_idx] += vap_oil_b;
1277  }
1278 
1279  for (int p = 0; p < this->number_of_phases_; ++p) {
1280  // TODO: double check the indices here
1281  this->ipr_a_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_a_perf[p];
1282  this->ipr_b_[this->ebosCompIdxToFlowCompIdx(p)] += ipr_b_perf[p];
1283  }
1284  }
1285  }
1286  }
1287 
1288  template<typename TypeTag>
1289  void
1290  MultisegmentWell<TypeTag>::
1291  checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& /*well_state*/, DeferredLogger& deferred_logger)
1292  {
1293  const auto& summaryState = ebos_simulator.vanguard().summaryState();
1294  const auto obtain_bhp = computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger);
1295 
1296  if (obtain_bhp) {
1297  this->operability_status_.can_obtain_bhp_with_thp_limit = true;
1298 
1299  const double bhp_limit = Base::mostStrictBhpFromBhpLimits(summaryState);
1300  this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1301 
1302  const double thp_limit = this->getTHPConstraint(summaryState);
1303  if (*obtain_bhp < thp_limit) {
1304  const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1305  + " bars is SMALLER than thp limit "
1306  + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1307  + " bars as a producer for well " + this->name();
1308  deferred_logger.debug(msg);
1309  }
1310  } else {
1311  // Shutting wells that can not find bhp value from thp
1312  // when under THP control
1313  this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1314  this->operability_status_.obey_bhp_limit_with_thp_limit = false;
1315  if (!this->wellIsStopped()) {
1316  const double thp_limit = this->getTHPConstraint(summaryState);
1317  deferred_logger.debug(" could not find bhp value at thp limit "
1318  + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1319  + " bar for well " + this->name() + ", the well might need to be closed ");
1320  }
1321  }
1322  }
1323 
1324 
1325 
1326 
1327 
1328  template<typename TypeTag>
1329  bool
1330  MultisegmentWell<TypeTag>::
1331  iterateWellEqWithControl(const Simulator& ebosSimulator,
1332  const double dt,
1333  const Well::InjectionControls& inj_controls,
1334  const Well::ProductionControls& prod_controls,
1335  WellState& well_state,
1336  const GroupState& group_state,
1337  DeferredLogger& deferred_logger)
1338  {
1339  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
1340 
1341  const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
1342  const WellState well_state0 = well_state;
1343  const std::vector<Scalar> residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger);
1344  std::vector<std::vector<Scalar> > residual_history;
1345  std::vector<double> measure_history;
1346  int it = 0;
1347  // relaxation factor
1348  double relaxation_factor = 1.;
1349  const double min_relaxation_factor = 0.6;
1350  bool converged = false;
1351  int stagnate_count = 0;
1352  bool relax_convergence = false;
1353  for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
1354 
1355  assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
1356 
1357  const BVectorWell dx_well = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
1358 
1359  if (it > this->param_.strict_inner_iter_wells_)
1360  relax_convergence = true;
1361 
1362  const auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
1363  if (report.converged()) {
1364  converged = true;
1365  break;
1366  }
1367 
1368  residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger));
1369  measure_history.push_back(this->getResidualMeasureValue(well_state,
1370  residual_history[it],
1371  this->param_.tolerance_wells_,
1372  this->param_.tolerance_pressure_ms_wells_,
1373  deferred_logger) );
1374 
1375  bool is_oscillate = false;
1376  bool is_stagnate = false;
1377 
1378  this->detectOscillations(measure_history, it, is_oscillate, is_stagnate);
1379  // TODO: maybe we should have more sophiscated strategy to recover the relaxation factor,
1380  // for example, to recover it to be bigger
1381 
1382  if (is_oscillate || is_stagnate) {
1383  // HACK!
1384  std::ostringstream sstr;
1385  if (relaxation_factor == min_relaxation_factor) {
1386  // Still stagnating, terminate iterations if 5 iterations pass.
1387  ++stagnate_count;
1388  if (stagnate_count == 6) {
1389  sstr << " well " << this->name() << " observes severe stagnation and/or oscillation. We relax the tolerance and check for convergence. \n";
1390  const auto reportStag = getWellConvergence(well_state, Base::B_avg_, deferred_logger, true);
1391  if (reportStag.converged()) {
1392  converged = true;
1393  sstr << " well " << this->name() << " manages to get converged with relaxed tolerances in " << it << " inner iterations";
1394  deferred_logger.debug(sstr.str());
1395  return converged;
1396  }
1397  }
1398  }
1399 
1400  // a factor value to reduce the relaxation_factor
1401  const double reduction_mutliplier = 0.9;
1402  relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
1403 
1404  // debug output
1405  if (is_stagnate) {
1406  sstr << " well " << this->name() << " observes stagnation in inner iteration " << it << "\n";
1407 
1408  }
1409  if (is_oscillate) {
1410  sstr << " well " << this->name() << " observes oscillation in inner iteration " << it << "\n";
1411  }
1412  sstr << " relaxation_factor is " << relaxation_factor << " now\n";
1413  deferred_logger.debug(sstr.str());
1414  }
1415  updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
1416  initPrimaryVariablesEvaluation();
1417  }
1418 
1419  // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state
1420  if (converged) {
1421  std::ostringstream sstr;
1422  sstr << " Well " << this->name() << " converged in " << it << " inner iterations.";
1423  if (relax_convergence)
1424  sstr << " (A relaxed tolerance was used after "<< this->param_.strict_inner_iter_wells_ << " iterations)";
1425  deferred_logger.debug(sstr.str());
1426  } else {
1427  std::ostringstream sstr;
1428  sstr << " Well " << this->name() << " did not converge in " << it << " inner iterations.";
1429 #define EXTRA_DEBUG_MSW 0
1430 #if EXTRA_DEBUG_MSW
1431  sstr << "***** Outputting the residual history for well " << this->name() << " during inner iterations:";
1432  for (int i = 0; i < it; ++i) {
1433  const auto& residual = residual_history[i];
1434  sstr << " residual at " << i << "th iteration ";
1435  for (const auto& res : residual) {
1436  sstr << " " << res;
1437  }
1438  sstr << " " << measure_history[i] << " \n";
1439  }
1440 #endif
1441  deferred_logger.debug(sstr.str());
1442  }
1443 
1444  return converged;
1445  }
1446 
1447 
1448 
1449 
1450 
1451  template<typename TypeTag>
1452  void
1453  MultisegmentWell<TypeTag>::
1454  assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
1455  const double dt,
1456  const Well::InjectionControls& inj_controls,
1457  const Well::ProductionControls& prod_controls,
1458  WellState& well_state,
1459  const GroupState& group_state,
1460  DeferredLogger& deferred_logger)
1461  {
1462 
1463  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1464 
1465  // update the upwinding segments
1466  this->updateUpwindingSegments();
1467 
1468  // calculate the fluid properties needed.
1469  computeSegmentFluidProperties(ebosSimulator);
1470 
1471  // clear all entries
1472  this->duneB_ = 0.0;
1473  this->duneC_ = 0.0;
1474 
1475  this->duneD_ = 0.0;
1476  this->resWell_ = 0.0;
1477 
1478  this->duneDSolver_.reset();
1479 
1480  auto& ws = well_state.well(this->index_of_well_);
1481  ws.dissolved_gas_rate = 0;
1482  ws.vaporized_oil_rate = 0;
1483 
1484  // for the black oil cases, there will be four equations,
1485  // the first three of them are the mass balance equations, the last one is the pressure equations.
1486  //
1487  // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same.
1488 
1489  const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1490 
1491  const int nseg = this->numberOfSegments();
1492 
1493  for (int seg = 0; seg < nseg; ++seg) {
1494  // calculating the accumulation term
1495  // TODO: without considering the efficiencty factor for now
1496  {
1497  const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
1498 
1499  // Add a regularization_factor to increase the accumulation term
1500  // This will make the system less stiff and help convergence for
1501  // difficult cases
1502  const Scalar regularization_factor = this->param_.regularization_factor_ms_wells_;
1503  // for each component
1504  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1505  const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->surfaceVolumeFraction(seg, comp_idx)
1506  - segment_fluid_initial_[seg][comp_idx]) / dt;
1507 
1508  this->resWell_[seg][comp_idx] += accumulation_term.value();
1509  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1510  this->duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + Indices::numEq);
1511  }
1512  }
1513  }
1514  // considering the contributions due to flowing out from the segment
1515  {
1516  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1517  const EvalWell segment_rate = this->getSegmentRateUpwinding(seg, comp_idx) * this->well_efficiency_factor_;
1518 
1519  const int seg_upwind = this->upwinding_segments_[seg];
1520  // segment_rate contains the derivatives with respect to GTotal in seg,
1521  // and WFrac and GFrac in seg_upwind
1522  this->resWell_[seg][comp_idx] -= segment_rate.value();
1523  this->duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + Indices::numEq);
1524  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1525  this->duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + Indices::numEq);
1526  }
1527  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1528  this->duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + Indices::numEq);
1529  }
1530  // pressure derivative should be zero
1531  }
1532  }
1533 
1534  // considering the contributions from the inlet segments
1535  {
1536  for (const int inlet : this->segment_inlets_[seg]) {
1537  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1538  const EvalWell inlet_rate = this->getSegmentRateUpwinding(inlet, comp_idx) * this->well_efficiency_factor_;
1539 
1540  const int inlet_upwind = this->upwinding_segments_[inlet];
1541  // inlet_rate contains the derivatives with respect to GTotal in inlet,
1542  // and WFrac and GFrac in inlet_upwind
1543  this->resWell_[seg][comp_idx] += inlet_rate.value();
1544  this->duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + Indices::numEq);
1545  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
1546  this->duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + Indices::numEq);
1547  }
1548  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1549  this->duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + Indices::numEq);
1550  }
1551  // pressure derivative should be zero
1552  }
1553  }
1554  }
1555 
1556  // calculating the perforation rate for each perforation that belongs to this segment
1557  const EvalWell seg_pressure = this->getSegmentPressure(seg);
1558  auto& perf_data = ws.perf_data;
1559  auto& perf_rates = perf_data.phase_rates;
1560  auto& perf_press_state = perf_data.pressure;
1561  for (const int perf : this->segment_perforations_[seg]) {
1562  const int cell_idx = this->well_cells_[perf];
1563  const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1564  std::vector<EvalWell> mob(this->num_components_, 0.0);
1565  getMobilityEval(ebosSimulator, perf, mob);
1566  const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1567  const double Tw = this->well_index_[perf] * trans_mult;
1568  std::vector<EvalWell> cq_s(this->num_components_, 0.0);
1569  EvalWell perf_press;
1570  double perf_dis_gas_rate = 0.;
1571  double perf_vap_oil_rate = 0.;
1572  computePerfRateEval(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
1573 
1574  // updating the solution gas rate and solution oil rate
1575  if (this->isProducer()) {
1576  ws.dissolved_gas_rate += perf_dis_gas_rate;
1577  ws.vaporized_oil_rate += perf_vap_oil_rate;
1578  }
1579 
1580  // store the perf pressure and rates
1581  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1582  perf_rates[perf*this->number_of_phases_ + this->ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
1583  }
1584  perf_press_state[perf] = perf_press.value();
1585 
1586  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1587  // the cq_s entering mass balance equations need to consider the efficiency factors.
1588  const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
1589 
1590  this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective);
1591 
1592  // subtract sum of phase fluxes in the well equations.
1593  this->resWell_[seg][comp_idx] += cq_s_effective.value();
1594 
1595  // assemble the jacobians
1596  for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
1597 
1598  // also need to consider the efficiency factor when manipulating the jacobians.
1599  this->duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + Indices::numEq); // intput in transformed matrix
1600 
1601  // the index name for the D should be eq_idx / pv_idx
1602  this->duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + Indices::numEq);
1603  }
1604 
1605  for (int pv_idx = 0; pv_idx < Indices::numEq; ++pv_idx) {
1606  // also need to consider the efficiency factor when manipulating the jacobians.
1607  this->duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx);
1608  }
1609  }
1610  }
1611 
1612  // the fourth dequation, the pressure drop equation
1613  if (seg == 0) { // top segment, pressure equation is the control equation
1614  const auto& summaryState = ebosSimulator.vanguard().summaryState();
1615  const Schedule& schedule = ebosSimulator.vanguard().schedule();
1616  this->assembleControlEq(well_state,
1617  group_state,
1618  schedule,
1619  summaryState,
1620  inj_controls,
1621  prod_controls,
1622  getRefDensity(),
1623  deferred_logger);
1624  } else {
1625  const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
1626  this->assemblePressureEq(seg, unit_system, well_state, deferred_logger);
1627  }
1628  }
1629  }
1630 
1631 
1632 
1633 
1634  template<typename TypeTag>
1635  bool
1636  MultisegmentWell<TypeTag>::
1637  openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
1638  {
1639  return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1640  }
1641 
1642 
1643  template<typename TypeTag>
1644  bool
1645  MultisegmentWell<TypeTag>::
1646  allDrawDownWrongDirection(const Simulator& ebos_simulator) const
1647  {
1648  bool all_drawdown_wrong_direction = true;
1649  const int nseg = this->numberOfSegments();
1650 
1651  for (int seg = 0; seg < nseg; ++seg) {
1652  const EvalWell segment_pressure = this->getSegmentPressure(seg);
1653  for (const int perf : this->segment_perforations_[seg]) {
1654 
1655  const int cell_idx = this->well_cells_[perf];
1656  const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1657  const auto& fs = intQuants.fluidState();
1658 
1659  // pressure difference between the segment and the perforation
1660  const EvalWell perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
1661  // pressure difference between the perforation and the grid cell
1662  const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
1663 
1664  const double pressure_cell = (fs.pressure(FluidSystem::oilPhaseIdx)).value();
1665  const double perf_press = pressure_cell - cell_perf_press_diff;
1666  // Pressure drawdown (also used to determine direction of flow)
1667  // TODO: not 100% sure about the sign of the seg_perf_press_diff
1668  const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
1669 
1670  // for now, if there is one perforation can produce/inject in the correct
1671  // direction, we consider this well can still produce/inject.
1672  // TODO: it can be more complicated than this to cause wrong-signed rates
1673  if ( (drawdown < 0. && this->isInjector()) ||
1674  (drawdown > 0. && this->isProducer()) ) {
1675  all_drawdown_wrong_direction = false;
1676  break;
1677  }
1678  }
1679  }
1680 
1681  return all_drawdown_wrong_direction;
1682  }
1683 
1684 
1685 
1686 
1687  template<typename TypeTag>
1688  void
1689  MultisegmentWell<TypeTag>::
1690  updateWaterThroughput(const double /*dt*/, WellState& /*well_state*/) const
1691  {
1692  }
1693 
1694 
1695 
1696 
1697 
1698  template<typename TypeTag>
1699  typename MultisegmentWell<TypeTag>::EvalWell
1700  MultisegmentWell<TypeTag>::
1701  getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const
1702  {
1703  EvalWell temperature;
1704  EvalWell saltConcentration;
1705  int pvt_region_index;
1706  {
1707  // using the pvt region of first perforated cell
1708  // TODO: it should be a member of the WellInterface, initialized properly
1709  const int cell_idx = this->well_cells_[0];
1710  const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1711  const auto& fs = intQuants.fluidState();
1712  temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value());
1713  saltConcentration = this->extendEval(fs.saltConcentration());
1714  pvt_region_index = fs.pvtRegionIndex();
1715  }
1716 
1717  return this->MSWEval::getSegmentSurfaceVolume(temperature,
1718  saltConcentration,
1719  pvt_region_index,
1720  seg_idx);
1721  }
1722 
1723 
1724 
1725 
1726 
1727  template<typename TypeTag>
1728  std::optional<double>
1729  MultisegmentWell<TypeTag>::
1730  computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
1731  const SummaryState& summary_state,
1732  DeferredLogger& deferred_logger) const
1733  {
1734  // Make the frates() function.
1735  auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1736  // Not solving the well equations here, which means we are
1737  // calculating at the current Fg/Fw values of the
1738  // well. This does not matter unless the well is
1739  // crossflowing, and then it is likely still a good
1740  // approximation.
1741  std::vector<double> rates(3);
1742  computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1743  return rates;
1744  };
1745 
1746  auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1747  computeBhpAtThpLimitProd(frates,
1748  summary_state,
1749  maxPerfPress(ebos_simulator),
1750  getRefDensity(),
1751  deferred_logger);
1752 
1753  if(bhpAtLimit)
1754  return bhpAtLimit;
1755 
1756  auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1757  // Solver the well iterations to see if we are
1758  // able to get a solution with an update
1759  // solution
1760  std::vector<double> rates(3);
1761  computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1762  return rates;
1763  };
1764 
1765  return this->MultisegmentWellGeneric<Scalar>::
1766  computeBhpAtThpLimitProd(fratesIter,
1767  summary_state,
1768  maxPerfPress(ebos_simulator),
1769  getRefDensity(),
1770  deferred_logger);
1771 
1772  }
1773 
1774 
1775 
1776 
1777  template<typename TypeTag>
1778  std::optional<double>
1779  MultisegmentWell<TypeTag>::
1780  computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
1781  const SummaryState& summary_state,
1782  DeferredLogger& deferred_logger) const
1783  {
1784  // Make the frates() function.
1785  auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1786  // Not solving the well equations here, which means we are
1787  // calculating at the current Fg/Fw values of the
1788  // well. This does not matter unless the well is
1789  // crossflowing, and then it is likely still a good
1790  // approximation.
1791  std::vector<double> rates(3);
1792  computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
1793  return rates;
1794  };
1795 
1796  auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
1797  computeBhpAtThpLimitInj(frates,
1798  summary_state,
1799  getRefDensity(),
1800  deferred_logger);
1801 
1802  if(bhpAtLimit)
1803  return bhpAtLimit;
1804 
1805  auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
1806  // Solver the well iterations to see if we are
1807  // able to get a solution with an update
1808  // solution
1809  std::vector<double> rates(3);
1810  computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
1811  return rates;
1812  };
1813 
1814  return this->MultisegmentWellGeneric<Scalar>::
1815  computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger);
1816  }
1817 
1818 
1819 
1820 
1821 
1822  template<typename TypeTag>
1823  double
1824  MultisegmentWell<TypeTag>::
1825  maxPerfPress(const Simulator& ebos_simulator) const
1826  {
1827  double max_pressure = 0.0;
1828  const int nseg = this->numberOfSegments();
1829  for (int seg = 0; seg < nseg; ++seg) {
1830  for (const int perf : this->segment_perforations_[seg]) {
1831  const int cell_idx = this->well_cells_[perf];
1832  const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1833  const auto& fs = int_quants.fluidState();
1834  double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value();
1835  max_pressure = std::max(max_pressure, pressure_cell);
1836  }
1837  }
1838  return max_pressure;
1839  }
1840 
1841 
1842 
1843 
1844 
1845  template<typename TypeTag>
1846  std::vector<double>
1848  computeCurrentWellRates(const Simulator& ebosSimulator,
1849  DeferredLogger& deferred_logger) const
1850  {
1851  // Calculate the rates that follow from the current primary variables.
1852  std::vector<Scalar> well_q_s(this->num_components_, 0.0);
1853  const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
1854  const int nseg = this->numberOfSegments();
1855  for (int seg = 0; seg < nseg; ++seg) {
1856  // calculating the perforation rate for each perforation that belongs to this segment
1857  const Scalar seg_pressure = getValue(this->getSegmentPressure(seg));
1858  for (const int perf : this->segment_perforations_[seg]) {
1859  const int cell_idx = this->well_cells_[perf];
1860  const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1861  std::vector<Scalar> mob(this->num_components_, 0.0);
1862  getMobilityScalar(ebosSimulator, perf, mob);
1863  const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
1864  const double Tw = this->well_index_[perf] * trans_mult;
1865  std::vector<Scalar> cq_s(this->num_components_, 0.0);
1866  computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, deferred_logger);
1867  for (int comp = 0; comp < this->num_components_; ++comp) {
1868  well_q_s[comp] += cq_s[comp];
1869  }
1870  }
1871  }
1872  return well_q_s;
1873  }
1874 
1875 
1876 
1877 
1878 
1879  template<typename TypeTag>
1880  void
1882  computeConnLevelProdInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
1883  const std::function<double(const double)>& connPICalc,
1884  const std::vector<Scalar>& mobility,
1885  double* connPI) const
1886  {
1887  const auto& pu = this->phaseUsage();
1888  const int np = this->number_of_phases_;
1889  for (int p = 0; p < np; ++p) {
1890  // Note: E100's notion of PI value phase mobility includes
1891  // the reciprocal FVF.
1892  const auto connMob =
1893  mobility[ this->flowPhaseToEbosCompIdx(p) ]
1894  * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
1895 
1896  connPI[p] = connPICalc(connMob);
1897  }
1898 
1899  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
1900  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
1901  {
1902  const auto io = pu.phase_pos[Oil];
1903  const auto ig = pu.phase_pos[Gas];
1904 
1905  const auto vapoil = connPI[ig] * fs.Rv().value();
1906  const auto disgas = connPI[io] * fs.Rs().value();
1907 
1908  connPI[io] += vapoil;
1909  connPI[ig] += disgas;
1910  }
1911  }
1912 
1913 
1914 
1915 
1916 
1917  template<typename TypeTag>
1918  void
1919  MultisegmentWell<TypeTag>::
1920  computeConnLevelInjInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
1921  const Phase preferred_phase,
1922  const std::function<double(const double)>& connIICalc,
1923  const std::vector<Scalar>& mobility,
1924  double* connII,
1925  DeferredLogger& deferred_logger) const
1926  {
1927  // Assumes single phase injection
1928  const auto& pu = this->phaseUsage();
1929 
1930  auto phase_pos = 0;
1931  if (preferred_phase == Phase::GAS) {
1932  phase_pos = pu.phase_pos[Gas];
1933  }
1934  else if (preferred_phase == Phase::OIL) {
1935  phase_pos = pu.phase_pos[Oil];
1936  }
1937  else if (preferred_phase == Phase::WATER) {
1938  phase_pos = pu.phase_pos[Water];
1939  }
1940  else {
1941  OPM_DEFLOG_THROW(NotImplemented,
1942  "Unsupported Injector Type ("
1943  << static_cast<int>(preferred_phase)
1944  << ") for well " << this->name()
1945  << " during connection I.I. calculation",
1946  deferred_logger);
1947  }
1948 
1949  const Scalar mt = std::accumulate(mobility.begin(), mobility.end(), 0.0);
1950  connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
1951  }
1952 
1953 } // namespace Opm
Represents the convergence status of the whole simulator, to make it possible to query and store the ...
Definition: ConvergenceReport.hpp:36
Definition: DeferredLogger.hpp:57
Definition: GroupState.hpp:34
Definition: MultisegmentWell.hpp:39
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:26
PhaseUsage phaseUsage(const Phases &phases)
Determine the active phases.
Definition: phaseUsageFromDeck.cpp:33