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