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