My Project
StandardWell_impl.hpp
1 /*
2  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
3  Copyright 2017 Statoil ASA.
4  Copyright 2016 - 2017 IRIS AS.
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #include <opm/common/utility/numeric/RootFinders.hpp>
23 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
24 #include <opm/simulators/linalg/SmallDenseMatrixUtils.hpp>
25 #include <opm/simulators/wells/VFPHelpers.hpp>
26 
27 #include <algorithm>
28 #include <functional>
29 #include <numeric>
30 
31 namespace Opm
32 {
33 
34  template<typename TypeTag>
35  StandardWell<TypeTag>::
36  StandardWell(const Well& well,
37  const ParallelWellInfo& pw_info,
38  const int time_step,
39  const ModelParameters& param,
40  const RateConverterType& rate_converter,
41  const int pvtRegionIdx,
42  const int num_components,
43  const int num_phases,
44  const int index_of_well,
45  const std::vector<PerforationData>& perf_data)
46  : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
47  , StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices,Scalar>&>(*this))
48  , regularize_(false)
49  {
50  assert(this->num_components_ == numWellConservationEq);
51  }
52 
53 
54 
55 
56 
57  template<typename TypeTag>
58  void
59  StandardWell<TypeTag>::
60  init(const PhaseUsage* phase_usage_arg,
61  const std::vector<double>& depth_arg,
62  const double gravity_arg,
63  const int num_cells,
64  const std::vector< Scalar >& B_avg,
65  const bool changed_to_open_this_step)
66  {
67  Base::init(phase_usage_arg, depth_arg, gravity_arg, num_cells, B_avg, changed_to_open_this_step);
68  this->StdWellEval::init(this->perf_depth_, depth_arg, num_cells, Base::has_polymermw);
69  }
70 
71 
72 
73 
74 
75  template<typename TypeTag>
76  void StandardWell<TypeTag>::
77  initPrimaryVariablesEvaluation() const
78  {
79  this->StdWellEval::initPrimaryVariablesEvaluation();
80  }
81 
82 
83 
84 
85 
86  template<typename TypeTag>
87  void
88  StandardWell<TypeTag>::
89  computePerfRateEval(const IntensiveQuantities& intQuants,
90  const std::vector<EvalWell>& mob,
91  const EvalWell& bhp,
92  const double Tw,
93  const int perf,
94  const bool allow_cf,
95  std::vector<EvalWell>& cq_s,
96  double& perf_dis_gas_rate,
97  double& perf_vap_oil_rate,
98  double& perf_vap_wat_rate,
99  DeferredLogger& deferred_logger) const
100  {
101  const auto& fs = intQuants.fluidState();
102  const EvalWell pressure = this->extendEval(this->getPerfCellPressure(fs));
103  const EvalWell rs = this->extendEval(fs.Rs());
104  const EvalWell rv = this->extendEval(fs.Rv());
105  const EvalWell rvw = this->extendEval(fs.Rvw());
106 
107  std::vector<EvalWell> b_perfcells_dense(this->num_components_, EvalWell{this->numWellEq_ + Indices::numEq, 0.0});
108  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
109  if (!FluidSystem::phaseIsActive(phaseIdx)) {
110  continue;
111  }
112  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
113  b_perfcells_dense[compIdx] = this->extendEval(fs.invB(phaseIdx));
114  }
115  if constexpr (has_solvent) {
116  b_perfcells_dense[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventInverseFormationVolumeFactor());
117  }
118 
119  if constexpr (has_zFraction) {
120  if (this->isInjector()) {
121  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
122  b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
123  b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
124  }
125  }
126 
127  EvalWell skin_pressure = EvalWell{this->numWellEq_ + Indices::numEq, 0.0};
128  if (has_polymermw) {
129  if (this->isInjector()) {
130  const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
131  skin_pressure = this->primary_variables_evaluation_[pskin_index];
132  }
133  }
134 
135  // surface volume fraction of fluids within wellbore
136  std::vector<EvalWell> cmix_s(this->numComponents(), EvalWell{this->numWellEq_ + Indices::numEq});
137  for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
138  cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
139  }
140 
141  computePerfRate(mob,
142  pressure,
143  bhp,
144  rs,
145  rv,
146  rvw,
147  b_perfcells_dense,
148  Tw,
149  perf,
150  allow_cf,
151  skin_pressure,
152  cmix_s,
153  cq_s,
154  perf_dis_gas_rate,
155  perf_vap_oil_rate,
156  perf_vap_wat_rate,
157  deferred_logger);
158  }
159 
160  template<typename TypeTag>
161  void
162  StandardWell<TypeTag>::
163  computePerfRateScalar(const IntensiveQuantities& intQuants,
164  const std::vector<Scalar>& mob,
165  const Scalar& bhp,
166  const double Tw,
167  const int perf,
168  const bool allow_cf,
169  std::vector<Scalar>& cq_s,
170  DeferredLogger& deferred_logger) const
171  {
172  const auto& fs = intQuants.fluidState();
173  const Scalar pressure = this->getPerfCellPressure(fs).value();
174  const Scalar rs = fs.Rs().value();
175  const Scalar rv = fs.Rv().value();
176  const Scalar rvw = fs.Rvw().value();
177  std::vector<Scalar> b_perfcells_dense(this->num_components_, 0.0);
178  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
179  if (!FluidSystem::phaseIsActive(phaseIdx)) {
180  continue;
181  }
182  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
183  b_perfcells_dense[compIdx] = fs.invB(phaseIdx).value();
184  }
185  if constexpr (has_solvent) {
186  b_perfcells_dense[Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
187  }
188 
189  if constexpr (has_zFraction) {
190  if (this->isInjector()) {
191  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
192  b_perfcells_dense[gasCompIdx] *= (1.0 - this->wsolvent());
193  b_perfcells_dense[gasCompIdx] += this->wsolvent()*intQuants.zPureInvFormationVolumeFactor().value();
194  }
195  }
196 
197  Scalar skin_pressure =0.0;
198  if (has_polymermw) {
199  if (this->isInjector()) {
200  const int pskin_index = Bhp + 1 + this->numPerfs() + perf;
201  skin_pressure = getValue(this->primary_variables_evaluation_[pskin_index]);
202  }
203  }
204 
205  Scalar perf_dis_gas_rate = 0.0;
206  Scalar perf_vap_oil_rate = 0.0;
207  Scalar perf_vap_wat_rate = 0.0;
208 
209  // surface volume fraction of fluids within wellbore
210  std::vector<Scalar> cmix_s(this->numComponents(), 0.0);
211  for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
212  cmix_s[componentIdx] = getValue(this->wellSurfaceVolumeFraction(componentIdx));
213  }
214 
215  computePerfRate(mob,
216  pressure,
217  bhp,
218  rs,
219  rv,
220  rvw,
221  b_perfcells_dense,
222  Tw,
223  perf,
224  allow_cf,
225  skin_pressure,
226  cmix_s,
227  cq_s,
228  perf_dis_gas_rate,
229  perf_vap_oil_rate,
230  perf_vap_wat_rate,
231  deferred_logger);
232  }
233 
234  template<typename TypeTag>
235  template<class Value>
236  void
237  StandardWell<TypeTag>::
238  computePerfRate(const std::vector<Value>& mob,
239  const Value& pressure,
240  const Value& bhp,
241  const Value& rs,
242  const Value& rv,
243  const Value& rvw,
244  std::vector<Value>& b_perfcells_dense,
245  const double Tw,
246  const int perf,
247  const bool allow_cf,
248  const Value& skin_pressure,
249  const std::vector<Value>& cmix_s,
250  std::vector<Value>& cq_s,
251  double& perf_dis_gas_rate,
252  double& perf_vap_oil_rate,
253  double& perf_vap_wat_rate,
254  DeferredLogger& deferred_logger) const
255  {
256  // Pressure drawdown (also used to determine direction of flow)
257  const Value well_pressure = bhp + this->perf_pressure_diffs_[perf];
258  Value drawdown = pressure - well_pressure;
259  if (this->isInjector()) {
260  drawdown += skin_pressure;
261  }
262 
263  // producing perforations
264  if ( drawdown > 0 ) {
265  //Do nothing if crossflow is not allowed
266  if (!allow_cf && this->isInjector()) {
267  return;
268  }
269 
270  // compute component volumetric rates at standard conditions
271  for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
272  const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
273  cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
274  }
275 
276  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
277  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
278  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
279  const Value cq_sOil = cq_s[oilCompIdx];
280  const Value cq_sGas = cq_s[gasCompIdx];
281  const Value dis_gas = rs * cq_sOil;
282  const Value vap_oil = rv * cq_sGas;
283 
284  cq_s[gasCompIdx] += dis_gas;
285  cq_s[oilCompIdx] += vap_oil;
286 
287  // recording the perforation solution gas rate and solution oil rates
288  if (this->isProducer()) {
289  perf_dis_gas_rate = getValue(dis_gas);
290  perf_vap_oil_rate = getValue(vap_oil);
291  }
292  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
293  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
294  const Value vap_wat = rvw * cq_sGas;
295  cq_s[waterCompIdx] += vap_wat;
296  if (this->isProducer())
297  perf_vap_wat_rate = getValue(vap_wat);
298  }
299  }
300 
301  } else {
302  //Do nothing if crossflow is not allowed
303  if (!allow_cf && this->isProducer()) {
304  return;
305  }
306 
307  // Using total mobilities
308  Value total_mob_dense = mob[0];
309  for (int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
310  total_mob_dense += mob[componentIdx];
311  }
312 
313  // injection perforations total volume rates
314  const Value cqt_i = - Tw * (total_mob_dense * drawdown);
315 
316  // compute volume ratio between connection at standard conditions
317  Value volumeRatio = bhp * 0.0; // initialize it with the correct type
318 ;
319  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
320  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
321  volumeRatio += cmix_s[waterCompIdx] / b_perfcells_dense[waterCompIdx];
322  }
323 
324  if constexpr (Indices::enableSolvent) {
325  volumeRatio += cmix_s[Indices::contiSolventEqIdx] / b_perfcells_dense[Indices::contiSolventEqIdx];
326  }
327 
328  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
329  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
330  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
331  // Incorporate RS/RV factors if both oil and gas active
332  const Value d = 1.0 - rv * rs;
333 
334  if (d <= 0.0) {
335  std::ostringstream sstr;
336  sstr << "Problematic d value " << d << " obtained for well " << this->name()
337  << " during computePerfRate calculations with rs " << rs
338  << ", rv " << rv << " and pressure " << pressure
339  << " obtaining d " << d
340  << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
341  << " for this connection.";
342  deferred_logger.debug(sstr.str());
343  }
344  const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
345  volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
346 
347  const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
348  volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
349  }
350  else {
351  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
352  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
353  volumeRatio += cmix_s[oilCompIdx] / b_perfcells_dense[oilCompIdx];
354  }
355  if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
356  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
357  volumeRatio += cmix_s[gasCompIdx] / b_perfcells_dense[gasCompIdx];
358  }
359  }
360 
361  // injecting connections total volumerates at standard conditions
362  Value cqt_is = cqt_i/volumeRatio;
363  for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
364  cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
365  }
366 
367  // calculating the perforation solution gas rate and solution oil rates
368  if (this->isProducer()) {
369  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
370  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
371  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
372  // TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
373  // s means standard condition, r means reservoir condition
374  // q_os = q_or * b_o + rv * q_gr * b_g
375  // q_gs = q_gr * b_g + rs * q_or * b_o
376  // d = 1.0 - rs * rv
377  // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
378  // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
379 
380  const double d = 1.0 - getValue(rv) * getValue(rs);
381 
382  if (d <= 0.0) {
383  std::ostringstream sstr;
384  sstr << "Problematic d value " << d << " obtained for well " << this->name()
385  << " during computePerfRate calculations with rs " << rs
386  << ", rv " << rv << " and pressure " << pressure
387  << " obtaining d " << d
388  << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
389  << " for this connection.";
390  deferred_logger.debug(sstr.str());
391  } else {
392  // vaporized oil into gas
393  // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
394  perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
395  // dissolved of gas in oil
396  // rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
397  perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
398  }
399  if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
400  // q_ws = q_wr * b_w + rvw * q_gr * b_g
401  // q_wr = 1 / b_w * (q_ws - rvw * q_gr * b_g) = 1 / b_w * (q_ws - rvw * 1 / d (q_gs - rs * q_os))
402  // vaporized water in gas
403  // rvw * q_gr * b_g = q_ws -q_wr *b_w = rvw * (q_gs -rs *q_os) / d
404  perf_vap_wat_rate = getValue(rvw) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
405  }
406  }
407  else if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
408  //no oil
409  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
410  perf_vap_wat_rate = getValue(rvw) * getValue(cq_s[gasCompIdx]);
411  }
412  }
413  }
414  }
415 
416 
417  template<typename TypeTag>
418  void
419  StandardWell<TypeTag>::
420  assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
421  const double dt,
422  const Well::InjectionControls& /*inj_controls*/,
423  const Well::ProductionControls& /*prod_controls*/,
424  WellState& well_state,
425  const GroupState& group_state,
426  DeferredLogger& deferred_logger)
427  {
428  // TODO: only_wells should be put back to save some computation
429  // for example, the matrices B C does not need to update if only_wells
430  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
431 
432  // clear all entries
433  this->duneB_ = 0.0;
434  this->duneC_ = 0.0;
435  this->duneD_ = 0.0;
436  this->resWell_ = 0.0;
437 
438  assembleWellEqWithoutIterationImpl(ebosSimulator, dt, well_state, group_state, deferred_logger);
439  }
440 
441 
442 
443 
444  template<typename TypeTag>
445  void
446  StandardWell<TypeTag>::
447  assembleWellEqWithoutIterationImpl(const Simulator& ebosSimulator,
448  const double dt,
449  WellState& well_state,
450  const GroupState& group_state,
451  DeferredLogger& deferred_logger)
452  {
453  // try to regularize equation if the well does not converge
454  const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
455  const double volume = 0.1 * unit::cubic(unit::feet) * regularization_factor;
456 
457  auto& ws = well_state.well(this->index_of_well_);
458 
459  ws.vaporized_oil_rate = 0;
460  ws.dissolved_gas_rate = 0;
461  ws.vaporized_wat_rate = 0;
462 
463  const int np = this->number_of_phases_;
464 
465  std::vector<RateVector> connectionRates = this->connectionRates_; // Copy to get right size.
466  auto& perf_data = ws.perf_data;
467  auto& perf_rates = perf_data.phase_rates;
468  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
469  // Calculate perforation quantities.
470  std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
471  EvalWell water_flux_s{this->numWellEq_ + Indices::numEq, 0.0};
472  EvalWell cq_s_zfrac_effective{this->numWellEq_ + Indices::numEq, 0.0};
473  calculateSinglePerf(ebosSimulator, perf, well_state, connectionRates, cq_s, water_flux_s, cq_s_zfrac_effective, deferred_logger);
474 
475  // Equation assembly for this perforation.
476  if constexpr (has_polymer && Base::has_polymermw) {
477  if (this->isInjector()) {
478  handleInjectivityEquations(ebosSimulator, well_state, perf, water_flux_s, deferred_logger);
479  }
480  }
481  const int cell_idx = this->well_cells_[perf];
482  for (int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
483  // the cq_s entering mass balance equations need to consider the efficiency factors.
484  const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
485 
486  connectionRates[perf][componentIdx] = Base::restrictEval(cq_s_effective);
487 
488  // subtract sum of phase fluxes in the well equations.
489  this->resWell_[0][componentIdx] += cq_s_effective.value();
490 
491  // assemble the jacobians
492  for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
493  // also need to consider the efficiency factor when manipulating the jacobians.
494  this->duneC_[0][cell_idx][pvIdx][componentIdx] -= cq_s_effective.derivative(pvIdx+Indices::numEq); // intput in transformed matrix
495  this->duneD_[0][0][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx+Indices::numEq);
496  }
497 
498  for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
499  this->duneB_[0][cell_idx][componentIdx][pvIdx] += cq_s_effective.derivative(pvIdx);
500  }
501 
502  // Store the perforation phase flux for later usage.
503  if (has_solvent && componentIdx == Indices::contiSolventEqIdx) {
504  auto& perf_rate_solvent = perf_data.solvent_rates;
505  perf_rate_solvent[perf] = cq_s[componentIdx].value();
506  } else {
507  perf_rates[perf*np + this->ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
508  }
509  }
510 
511  if constexpr (has_zFraction) {
512  for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
513  this->duneC_[0][cell_idx][pvIdx][Indices::contiZfracEqIdx] -= cq_s_zfrac_effective.derivative(pvIdx+Indices::numEq);
514  }
515  }
516  }
517  // Update the connection
518  this->connectionRates_ = connectionRates;
519 
520  // Accumulate dissolved gas and vaporized oil flow rates across all
521  // ranks sharing this well (this->index_of_well_).
522  {
523  const auto& comm = this->parallel_well_info_.communication();
524  ws.dissolved_gas_rate = comm.sum(ws.dissolved_gas_rate);
525  ws.vaporized_oil_rate = comm.sum(ws.vaporized_oil_rate);
526  ws.vaporized_wat_rate = comm.sum(ws.vaporized_wat_rate);
527  }
528 
529  // accumulate resWell_ and duneD_ in parallel to get effects of all perforations (might be distributed)
530  wellhelpers::sumDistributedWellEntries(this->duneD_[0][0], this->resWell_[0],
531  this->parallel_well_info_.communication());
532  // add vol * dF/dt + Q to the well equations;
533  for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
534  // TODO: following the development in MSW, we need to convert the volume of the wellbore to be surface volume
535  // since all the rates are under surface condition
536  EvalWell resWell_loc(this->numWellEq_ + Indices::numEq, 0.0);
537  if (FluidSystem::numActivePhases() > 1) {
538  assert(dt > 0);
539  resWell_loc += (this->wellSurfaceVolumeFraction(componentIdx) - this->F0_[componentIdx]) * volume / dt;
540  }
541  resWell_loc -= this->getQs(componentIdx) * this->well_efficiency_factor_;
542  for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
543  this->duneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+Indices::numEq);
544  }
545  this->resWell_[0][componentIdx] += resWell_loc.value();
546  }
547 
548  const auto& summaryState = ebosSimulator.vanguard().summaryState();
549  const Schedule& schedule = ebosSimulator.vanguard().schedule();
550  this->assembleControlEq(well_state, group_state, schedule, summaryState, deferred_logger);
551 
552 
553  // do the local inversion of D.
554  try {
555  this->invDuneD_ = this->duneD_; // Not strictly need if not cpr with well contributions is used
556  detail::invertMatrix(this->invDuneD_[0][0]);
557  } catch (NumericalProblem&) {
558  // for singular matrices, use identity as the inverse
559  this->invDuneD_[0][0] = 0.0;
560  for (size_t i = 0; i < this->invDuneD_[0][0].rows(); ++i) {
561  this->invDuneD_[0][0][i][i] = 1.0;
562  }
563  } catch( ... ) {
564  OPM_DEFLOG_THROW(NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
565  }
566  }
567 
568 
569 
570 
571  template<typename TypeTag>
572  void
573  StandardWell<TypeTag>::
574  calculateSinglePerf(const Simulator& ebosSimulator,
575  const int perf,
576  WellState& well_state,
577  std::vector<RateVector>& connectionRates,
578  std::vector<EvalWell>& cq_s,
579  EvalWell& water_flux_s,
580  EvalWell& cq_s_zfrac_effective,
581  DeferredLogger& deferred_logger) const
582  {
583  const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
584  const EvalWell& bhp = this->getBhp();
585  const int cell_idx = this->well_cells_[perf];
586  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
587  std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
588  getMobilityEval(ebosSimulator, perf, mob, deferred_logger);
589 
590  double perf_dis_gas_rate = 0.;
591  double perf_vap_oil_rate = 0.;
592  double perf_vap_wat_rate = 0.;
593  double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
594  const double Tw = this->well_index_[perf] * trans_mult;
595  computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
596  cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
597 
598  auto& ws = well_state.well(this->index_of_well_);
599  auto& perf_data = ws.perf_data;
600  if constexpr (has_polymer && Base::has_polymermw) {
601  if (this->isInjector()) {
602  // Store the original water flux computed from the reservoir quantities.
603  // It will be required to assemble the injectivity equations.
604  const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
605  water_flux_s = cq_s[water_comp_idx];
606  // Modify the water flux for the rest of this function to depend directly on the
607  // local water velocity primary variable.
608  handleInjectivityRate(ebosSimulator, perf, cq_s);
609  }
610  }
611 
612  // updating the solution gas rate and solution oil rate
613  if (this->isProducer()) {
614  ws.dissolved_gas_rate += perf_dis_gas_rate;
615  ws.vaporized_oil_rate += perf_vap_oil_rate;
616  ws.vaporized_wat_rate += perf_vap_wat_rate;
617  }
618 
619  if constexpr (has_energy) {
620  connectionRates[perf][Indices::contiEnergyEqIdx] = 0.0;
621  }
622 
623  if constexpr (has_energy) {
624 
625  auto fs = intQuants.fluidState();
626  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
627  if (!FluidSystem::phaseIsActive(phaseIdx)) {
628  continue;
629  }
630 
631  // convert to reservoir conditions
632  EvalWell cq_r_thermal(this->numWellEq_ + Indices::numEq, 0.);
633  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
634  const bool both_oil_gas = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
635  if ( !both_oil_gas || FluidSystem::waterPhaseIdx == phaseIdx ) {
636  cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
637  } else {
638  // remove dissolved gas and vapporized oil
639  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
640  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
641  // q_os = q_or * b_o + rv * q_gr * b_g
642  // q_gs = q_gr * g_g + rs * q_or * b_o
643  // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
644  // d = 1.0 - rs * rv
645  const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
646  if (d <= 0.0) {
647  std::ostringstream sstr;
648  sstr << "Problematic d value " << d << " obtained for well " << this->name()
649  << " during calculateSinglePerf with rs " << fs.Rs()
650  << ", rv " << fs.Rv()
651  << " obtaining d " << d
652  << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
653  << " for this connection.";
654  deferred_logger.debug(sstr.str());
655  cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
656  } else {
657  if(FluidSystem::gasPhaseIdx == phaseIdx) {
658  cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
659  } else if(FluidSystem::oilPhaseIdx == phaseIdx) {
660  // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
661  cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
662  }
663  }
664  }
665 
666  // change temperature for injecting fluids
667  if (this->isInjector() && cq_s[activeCompIdx] > 0.0){
668  // only handles single phase injection now
669  assert(this->well_ecl_.injectorType() != InjectorType::MULTI);
670  fs.setTemperature(this->well_ecl_.temperature());
671  typedef typename std::decay<decltype(fs)>::type::Scalar FsScalar;
672  typename FluidSystem::template ParameterCache<FsScalar> paramCache;
673  const unsigned pvtRegionIdx = intQuants.pvtRegionIndex();
674  paramCache.setRegionIndex(pvtRegionIdx);
675  paramCache.setMaxOilSat(ebosSimulator.problem().maxOilSaturation(cell_idx));
676  paramCache.updatePhase(fs, phaseIdx);
677 
678  const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
679  fs.setDensity(phaseIdx, rho);
680  const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
681  fs.setEnthalpy(phaseIdx, h);
682  }
683  // compute the thermal flux
684  cq_r_thermal *= this->extendEval(fs.enthalpy(phaseIdx)) * this->extendEval(fs.density(phaseIdx));
685  connectionRates[perf][Indices::contiEnergyEqIdx] += Base::restrictEval(cq_r_thermal);
686  }
687  }
688 
689  if constexpr (has_polymer) {
690  // TODO: the application of well efficiency factor has not been tested with an example yet
691  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
692  EvalWell cq_s_poly = cq_s[waterCompIdx];
693  if (this->isInjector()) {
694  cq_s_poly *= this->wpolymer();
695  } else {
696  cq_s_poly *= this->extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
697  }
698  // Note. Efficiency factor is handled in the output layer
699  auto& perf_rate_polymer = perf_data.polymer_rates;
700  perf_rate_polymer[perf] = cq_s_poly.value();
701 
702  cq_s_poly *= this->well_efficiency_factor_;
703  connectionRates[perf][Indices::contiPolymerEqIdx] = Base::restrictEval(cq_s_poly);
704 
705  if constexpr (Base::has_polymermw) {
706  updateConnectionRatePolyMW(cq_s_poly, intQuants, well_state, perf, connectionRates, deferred_logger);
707  }
708  }
709 
710  if constexpr (has_foam) {
711  // TODO: the application of well efficiency factor has not been tested with an example yet
712  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
713  EvalWell cq_s_foam = cq_s[gasCompIdx] * this->well_efficiency_factor_;
714  if (this->isInjector()) {
715  cq_s_foam *= this->wfoam();
716  } else {
717  cq_s_foam *= this->extendEval(intQuants.foamConcentration());
718  }
719  connectionRates[perf][Indices::contiFoamEqIdx] = Base::restrictEval(cq_s_foam);
720  }
721 
722  if constexpr (has_zFraction) {
723  // TODO: the application of well efficiency factor has not been tested with an example yet
724  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
725  cq_s_zfrac_effective = cq_s[gasCompIdx];
726  if (this->isInjector()) {
727  cq_s_zfrac_effective *= this->wsolvent();
728  } else if (cq_s_zfrac_effective.value() != 0.0) {
729  const double dis_gas_frac = perf_dis_gas_rate / cq_s_zfrac_effective.value();
730  cq_s_zfrac_effective *= this->extendEval(dis_gas_frac*intQuants.xVolume() + (1.0-dis_gas_frac)*intQuants.yVolume());
731  }
732  auto& perf_rate_solvent = perf_data.solvent_rates;
733  perf_rate_solvent[perf] = cq_s_zfrac_effective.value();
734 
735  cq_s_zfrac_effective *= this->well_efficiency_factor_;
736  connectionRates[perf][Indices::contiZfracEqIdx] = Base::restrictEval(cq_s_zfrac_effective);
737  }
738 
739  if constexpr (has_brine) {
740  // TODO: the application of well efficiency factor has not been tested with an example yet
741  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
742  // Correction salt rate; evaporated water does not contain salt
743  EvalWell cq_s_sm = cq_s[waterCompIdx] - perf_vap_wat_rate;
744  if (this->isInjector()) {
745  cq_s_sm *= this->wsalt();
746  } else {
747  cq_s_sm *= this->extendEval(intQuants.fluidState().saltConcentration());
748  }
749  // Note. Efficiency factor is handled in the output layer
750  auto& perf_rate_brine = perf_data.brine_rates;
751  perf_rate_brine[perf] = cq_s_sm.value();
752 
753  cq_s_sm *= this->well_efficiency_factor_;
754  connectionRates[perf][Indices::contiBrineEqIdx] = Base::restrictEval(cq_s_sm);
755  }
756 
757  if constexpr (has_micp) {
758  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
759  EvalWell cq_s_microbe = cq_s[waterCompIdx];
760  if (this->isInjector()) {
761  cq_s_microbe *= this->wmicrobes();
762  } else {
763  cq_s_microbe *= this->extendEval(intQuants.microbialConcentration());
764  }
765  connectionRates[perf][Indices::contiMicrobialEqIdx] = Base::restrictEval(cq_s_microbe);
766  EvalWell cq_s_oxygen = cq_s[waterCompIdx];
767  if (this->isInjector()) {
768  cq_s_oxygen *= this->woxygen();
769  } else {
770  cq_s_oxygen *= this->extendEval(intQuants.oxygenConcentration());
771  }
772  connectionRates[perf][Indices::contiOxygenEqIdx] = Base::restrictEval(cq_s_oxygen);
773  EvalWell cq_s_urea = cq_s[waterCompIdx];
774  if (this->isInjector()) {
775  cq_s_urea *= this->wurea();
776  } else {
777  cq_s_urea *= this->extendEval(intQuants.ureaConcentration());
778  }
779  connectionRates[perf][Indices::contiUreaEqIdx] = Base::restrictEval(cq_s_urea);
780  }
781 
782  // Store the perforation pressure for later usage.
783  perf_data.pressure[perf] = ws.bhp + this->perf_pressure_diffs_[perf];
784  }
785 
786 
787 
788 
789  template<typename TypeTag>
790  void
791  StandardWell<TypeTag>::
792  getMobilityEval(const Simulator& ebosSimulator,
793  const int perf,
794  std::vector<EvalWell>& mob,
795  DeferredLogger& deferred_logger) const
796  {
797  const int cell_idx = this->well_cells_[perf];
798  assert (int(mob.size()) == this->num_components_);
799  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
800  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
801 
802  // either use mobility of the perforation cell or calcualte its own
803  // based on passing the saturation table index
804  const int satid = this->saturation_table_number_[perf] - 1;
805  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
806  if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
807 
808  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
809  if (!FluidSystem::phaseIsActive(phaseIdx)) {
810  continue;
811  }
812 
813  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
814  mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx));
815  }
816  if (has_solvent) {
817  mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility());
818  }
819  } else {
820 
821  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
822  std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
823  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
824 
825  // reset the satnumvalue back to original
826  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
827 
828  // compute the mobility
829  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
830  if (!FluidSystem::phaseIsActive(phaseIdx)) {
831  continue;
832  }
833 
834  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
835  mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx));
836  }
837 
838  // this may not work if viscosity and relperms has been modified?
839  if constexpr (has_solvent) {
840  OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger);
841  }
842  }
843 
844  // modify the water mobility if polymer is present
845  if constexpr (has_polymer) {
846  if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
847  OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger);
848  }
849 
850  // for the cases related to polymer molecular weight, we assume fully mixing
851  // as a result, the polymer and water share the same viscosity
852  if constexpr (!Base::has_polymermw) {
853  updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
854  }
855  }
856  }
857 
858  template<typename TypeTag>
859  void
860  StandardWell<TypeTag>::
861  getMobilityScalar(const Simulator& ebosSimulator,
862  const int perf,
863  std::vector<Scalar>& mob,
864  DeferredLogger& deferred_logger) const
865  {
866  const int cell_idx = this->well_cells_[perf];
867  assert (int(mob.size()) == this->num_components_);
868  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
869  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
870 
871  // either use mobility of the perforation cell or calcualte its own
872  // based on passing the saturation table index
873  const int satid = this->saturation_table_number_[perf] - 1;
874  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
875  if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
876 
877  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
878  if (!FluidSystem::phaseIsActive(phaseIdx)) {
879  continue;
880  }
881 
882  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
883  mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx));
884  }
885  if (has_solvent) {
886  mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility());
887  }
888  } else {
889 
890  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
891  std::array<Eval,3> relativePerms = { 0.0, 0.0, 0.0 };
892  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
893 
894  // reset the satnumvalue back to original
895  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
896 
897  // compute the mobility
898  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
899  if (!FluidSystem::phaseIsActive(phaseIdx)) {
900  continue;
901  }
902 
903  const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
904  mob[activeCompIdx] = getValue(relativePerms[phaseIdx]) / getValue(intQuants.fluidState().viscosity(phaseIdx));
905  }
906 
907  // this may not work if viscosity and relperms has been modified?
908  if constexpr (has_solvent) {
909  OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger);
910  }
911  }
912 
913  // modify the water mobility if polymer is present
914  if constexpr (has_polymer) {
915  if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
916  OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger);
917  }
918 
919  // for the cases related to polymer molecular weight, we assume fully mixing
920  // as a result, the polymer and water share the same viscosity
921  if constexpr (!Base::has_polymermw) {
922  std::vector<EvalWell> mob_eval(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
923  updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger);
924  for (size_t i = 0; i < mob.size(); ++i) {
925  mob[i] = getValue(mob_eval[i]);
926  }
927  }
928  }
929  }
930 
931 
932 
933  template<typename TypeTag>
934  void
935  StandardWell<TypeTag>::
936  updateWellState(const BVectorWell& dwells,
937  WellState& well_state,
938  DeferredLogger& deferred_logger) const
939  {
940  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
941 
942  updatePrimaryVariablesNewton(dwells, well_state, deferred_logger);
943 
944  updateWellStateFromPrimaryVariables(well_state, deferred_logger);
945  Base::calculateReservoirRates(well_state.well(this->index_of_well_));
946  }
947 
948 
949 
950 
951 
952  template<typename TypeTag>
953  void
954  StandardWell<TypeTag>::
955  updatePrimaryVariablesNewton(const BVectorWell& dwells,
956  const WellState& /* well_state */,
957  DeferredLogger& deferred_logger) const
958  {
959  const double dFLimit = this->param_.dwell_fraction_max_;
960  const double dBHPLimit = this->param_.dbhp_max_rel_;
961  this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit);
962 
963  updateExtraPrimaryVariables(dwells);
964 
965  for (double v : this->primary_variables_) {
966  if(!isfinite(v))
967  OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after newton update well: " << this->name(), deferred_logger);
968  }
969 
970  }
971 
972 
973 
974 
975 
976  template<typename TypeTag>
977  void
978  StandardWell<TypeTag>::
979  updateExtraPrimaryVariables(const BVectorWell& dwells) const
980  {
981  // for the water velocity and skin pressure
982  if constexpr (Base::has_polymermw) {
983  this->updatePrimaryVariablesPolyMW(dwells);
984  }
985  }
986 
987 
988 
989 
990 
991  template<typename TypeTag>
992  void
993  StandardWell<TypeTag>::
994  updateWellStateFromPrimaryVariables(WellState& well_state, DeferredLogger& deferred_logger) const
995  {
996  this->StdWellEval::updateWellStateFromPrimaryVariables(well_state, deferred_logger);
997 
998  // other primary variables related to polymer injectivity study
999  if constexpr (Base::has_polymermw) {
1000  this->updateWellStateFromPrimaryVariablesPolyMW(well_state);
1001  }
1002  }
1003 
1004 
1005 
1006 
1007 
1008  template<typename TypeTag>
1009  void
1010  StandardWell<TypeTag>::
1011  updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const
1012  {
1013  // TODO: not handling solvent related here for now
1014 
1015  // initialize all the values to be zero to begin with
1016  std::fill(this->ipr_a_.begin(), this->ipr_a_.end(), 0.);
1017  std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
1018 
1019  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
1020  std::vector<Scalar> mob(this->num_components_, 0.0);
1021  getMobilityScalar(ebos_simulator, perf, mob, deferred_logger);
1022 
1023  const int cell_idx = this->well_cells_[perf];
1024  const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1025  const auto& fs = int_quantities.fluidState();
1026  // the pressure of the reservoir grid block the well connection is in
1027  double p_r = this->getPerfCellPressure(fs).value();
1028 
1029  // calculating the b for the connection
1030  std::vector<double> b_perf(this->num_components_);
1031  for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
1032  if (!FluidSystem::phaseIsActive(phase)) {
1033  continue;
1034  }
1035  const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
1036  b_perf[comp_idx] = fs.invB(phase).value();
1037  }
1038  if constexpr (has_solvent) {
1039  b_perf[Indices::contiSolventEqIdx] = int_quantities.solventInverseFormationVolumeFactor().value();
1040  }
1041 
1042  // the pressure difference between the connection and BHP
1043  const double h_perf = this->perf_pressure_diffs_[perf];
1044  const double pressure_diff = p_r - h_perf;
1045 
1046  // Let us add a check, since the pressure is calculated based on zero value BHP
1047  // it should not be negative anyway. If it is negative, we might need to re-formulate
1048  // to taking into consideration the crossflow here.
1049  if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
1050  deferred_logger.debug("CROSSFLOW_IPR",
1051  "cross flow found when updateIPR for well " + name()
1052  + " . The connection is ignored in IPR calculations");
1053  // we ignore these connections for now
1054  continue;
1055  }
1056 
1057  // the well index associated with the connection
1058  const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
1059 
1060  std::vector<double> ipr_a_perf(this->ipr_a_.size());
1061  std::vector<double> ipr_b_perf(this->ipr_b_.size());
1062  for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
1063  const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
1064  ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
1065  ipr_b_perf[comp_idx] += tw_mob;
1066  }
1067 
1068  // we need to handle the rs and rv when both oil and gas are present
1069  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
1070  const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1071  const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1072  const double rs = (fs.Rs()).value();
1073  const double rv = (fs.Rv()).value();
1074 
1075  const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
1076  const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
1077 
1078  ipr_a_perf[gas_comp_idx] += dis_gas_a;
1079  ipr_a_perf[oil_comp_idx] += vap_oil_a;
1080 
1081  const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
1082  const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
1083 
1084  ipr_b_perf[gas_comp_idx] += dis_gas_b;
1085  ipr_b_perf[oil_comp_idx] += vap_oil_b;
1086  }
1087 
1088  for (size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
1089  this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
1090  this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
1091  }
1092  }
1093  this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size());
1094  this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size());
1095  }
1096 
1097 
1098  template<typename TypeTag>
1099  void
1100  StandardWell<TypeTag>::
1101  checkOperabilityUnderBHPLimit(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger)
1102  {
1103  const auto& summaryState = ebos_simulator.vanguard().summaryState();
1104  const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1105  // Crude but works: default is one atmosphere.
1106  // TODO: a better way to detect whether the BHP is defaulted or not
1107  const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
1108  if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
1109  // if the BHP limit is not defaulted or the well does not have a THP limit
1110  // we need to check the BHP limit
1111  double total_ipr_mass_rate = 0.0;
1112  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1113  {
1114  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1115  continue;
1116  }
1117 
1118  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1119  const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
1120 
1121  const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
1122  total_ipr_mass_rate += ipr_rate * rho;
1123  }
1124  if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
1125  this->operability_status_.operable_under_only_bhp_limit = false;
1126  }
1127 
1128  // checking whether running under BHP limit will violate THP limit
1129  if (this->operability_status_.operable_under_only_bhp_limit && this->wellHasTHPConstraints(summaryState)) {
1130  // option 1: calculate well rates based on the BHP limit.
1131  // option 2: stick with the above IPR curve
1132  // we use IPR here
1133  std::vector<double> well_rates_bhp_limit;
1134  computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
1135 
1136  this->adaptRatesForVFP(well_rates_bhp_limit);
1137  const double thp = this->calculateThpFromBhp(well_state, well_rates_bhp_limit, bhp_limit, deferred_logger);
1138  const double thp_limit = this->getTHPConstraint(summaryState);
1139  if ( (this->isProducer() && thp < thp_limit) || (this->isInjector() && thp > thp_limit) ) {
1140  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1141  }
1142  }
1143  } else {
1144  // defaulted BHP and there is a THP constraint
1145  // default BHP limit is about 1 atm.
1146  // when applied the hydrostatic pressure correction dp,
1147  // most likely we get a negative value (bhp + dp)to search in the VFP table,
1148  // which is not desirable.
1149  // we assume we can operate under defaulted BHP limit and will violate the THP limit
1150  // when operating under defaulted BHP limit.
1151  this->operability_status_.operable_under_only_bhp_limit = true;
1152  this->operability_status_.obey_thp_limit_under_bhp_limit = false;
1153  }
1154  }
1155 
1156 
1157 
1158 
1159 
1160  template<typename TypeTag>
1161  void
1162  StandardWell<TypeTag>::
1163  checkOperabilityUnderTHPLimit(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger)
1164  {
1165  const auto& summaryState = ebos_simulator.vanguard().summaryState();
1166  const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger)
1167  : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
1168 
1169  if (obtain_bhp) {
1170  this->operability_status_.can_obtain_bhp_with_thp_limit = true;
1171 
1172  const double bhp_limit = this->mostStrictBhpFromBhpLimits(summaryState);
1173  this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
1174 
1175  const double thp_limit = this->getTHPConstraint(summaryState);
1176  if (this->isProducer() && *obtain_bhp < thp_limit) {
1177  const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1178  + " bars is SMALLER than thp limit "
1179  + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1180  + " bars as a producer for well " + name();
1181  deferred_logger.debug(msg);
1182  }
1183  else if (this->isInjector() && *obtain_bhp > thp_limit) {
1184  const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
1185  + " bars is LARGER than thp limit "
1186  + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1187  + " bars as a injector for well " + name();
1188  deferred_logger.debug(msg);
1189  }
1190  } else {
1191  this->operability_status_.can_obtain_bhp_with_thp_limit = false;
1192  this->operability_status_.obey_bhp_limit_with_thp_limit = false;
1193  if (!this->wellIsStopped()) {
1194  const double thp_limit = this->getTHPConstraint(summaryState);
1195  deferred_logger.debug(" could not find bhp value at thp limit "
1196  + std::to_string(unit::convert::to(thp_limit, unit::barsa))
1197  + " bar for well " + name() + ", the well might need to be closed ");
1198  }
1199  }
1200  }
1201 
1202 
1203 
1204 
1205 
1206  template<typename TypeTag>
1207  bool
1208  StandardWell<TypeTag>::
1209  allDrawDownWrongDirection(const Simulator& ebos_simulator) const
1210  {
1211  bool all_drawdown_wrong_direction = true;
1212 
1213  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
1214  const int cell_idx = this->well_cells_[perf];
1215  const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1216  const auto& fs = intQuants.fluidState();
1217 
1218  const double pressure = this->getPerfCellPressure(fs).value();
1219  const double bhp = this->getBhp().value();
1220 
1221  // Pressure drawdown (also used to determine direction of flow)
1222  const double well_pressure = bhp + this->perf_pressure_diffs_[perf];
1223  const double drawdown = pressure - well_pressure;
1224 
1225  // for now, if there is one perforation can produce/inject in the correct
1226  // direction, we consider this well can still produce/inject.
1227  // TODO: it can be more complicated than this to cause wrong-signed rates
1228  if ( (drawdown < 0. && this->isInjector()) ||
1229  (drawdown > 0. && this->isProducer()) ) {
1230  all_drawdown_wrong_direction = false;
1231  break;
1232  }
1233  }
1234 
1235  const auto& comm = this->parallel_well_info_.communication();
1236  if (comm.size() > 1)
1237  {
1238  all_drawdown_wrong_direction =
1239  (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1);
1240  }
1241 
1242  return all_drawdown_wrong_direction;
1243  }
1244 
1245 
1246 
1247 
1248  template<typename TypeTag>
1249  bool
1250  StandardWell<TypeTag>::
1251  canProduceInjectWithCurrentBhp(const Simulator& ebos_simulator,
1252  const WellState& well_state,
1253  DeferredLogger& deferred_logger)
1254  {
1255  const double bhp = well_state.well(this->index_of_well_).bhp;
1256  std::vector<double> well_rates;
1257  computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
1258 
1259  const double sign = (this->isProducer()) ? -1. : 1.;
1260  const double threshold = sign * std::numeric_limits<double>::min();
1261 
1262  bool can_produce_inject = false;
1263  for (const auto value : well_rates) {
1264  if (this->isProducer() && value < threshold) {
1265  can_produce_inject = true;
1266  break;
1267  } else if (this->isInjector() && value > threshold) {
1268  can_produce_inject = true;
1269  break;
1270  }
1271  }
1272 
1273  if (!can_produce_inject) {
1274  deferred_logger.debug(" well " + name() + " CANNOT produce or inejct ");
1275  }
1276 
1277  return can_produce_inject;
1278  }
1279 
1280 
1281 
1282 
1283 
1284  template<typename TypeTag>
1285  bool
1286  StandardWell<TypeTag>::
1287  openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const
1288  {
1289  return !this->getAllowCrossFlow() && allDrawDownWrongDirection(ebos_simulator);
1290  }
1291 
1292 
1293 
1294 
1295  template<typename TypeTag>
1296  void
1297  StandardWell<TypeTag>::
1298  computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
1299  const WellState& well_state,
1300  std::vector<double>& b_perf,
1301  std::vector<double>& rsmax_perf,
1302  std::vector<double>& rvmax_perf,
1303  std::vector<double>& rvwmax_perf,
1304  std::vector<double>& surf_dens_perf) const
1305  {
1306  const int nperf = this->number_of_perforations_;
1307  const PhaseUsage& pu = phaseUsage();
1308  b_perf.resize(nperf * this->num_components_);
1309  surf_dens_perf.resize(nperf * this->num_components_);
1310  const auto& ws = well_state.well(this->index_of_well_);
1311 
1312  const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
1313  const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
1314  const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
1315 
1316  //rs and rv are only used if both oil and gas is present
1317  if (oilPresent && gasPresent) {
1318  rsmax_perf.resize(nperf);
1319  rvmax_perf.resize(nperf);
1320  }
1321  //rvw is only used if both water and gas is present
1322  if (waterPresent && gasPresent) {
1323  rvwmax_perf.resize(nperf);
1324  }
1325 
1326  // Compute the average pressure in each well block
1327  const auto& perf_press = ws.perf_data.pressure;
1328  auto p_above = this->parallel_well_info_.communicateAboveValues(ws.bhp,
1329  perf_press.data(),
1330  nperf);
1331 
1332  for (int perf = 0; perf < nperf; ++perf) {
1333  const int cell_idx = this->well_cells_[perf];
1334  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1335  const auto& fs = intQuants.fluidState();
1336 
1337  const double p_avg = (perf_press[perf] + p_above[perf])/2;
1338  const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1339  const double saltConcentration = fs.saltConcentration().value();
1340 
1341  if (waterPresent) {
1342  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
1343  b_perf[ waterCompIdx + perf * this->num_components_] =
1344  FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, saltConcentration);
1345  }
1346 
1347  if (gasPresent) {
1348  const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
1349  const int gaspos = gasCompIdx + perf * this->num_components_;
1350 
1351  if (oilPresent && waterPresent) {
1352  const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
1353  const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers
1354  rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1355  rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1356  double rv = 0.0;
1357  double rvw = 0.0;
1358  if (oilrate > 0) {
1359  const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1360  if (gasrate > 0) {
1361  rv = oilrate / gasrate;
1362  }
1363  rv = std::min(rv, rvmax_perf[perf]);
1364  }
1365  if (waterrate > 0) {
1366  const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1367  if (gasrate > 0) {
1368  rvw = waterrate / gasrate;
1369  }
1370  rvw = std::min(rvw, rvwmax_perf[perf]);
1371  }
1372  if (rv > 0.0 || rvw > 0.0){
1373  b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, rvw);
1374  }
1375  else {
1376  b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1377  }
1378  } else if (oilPresent) {
1379  //no water
1380  const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
1381  rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1382  if (oilrate > 0) {
1383  const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1384  double rv = 0.0;
1385  if (gasrate > 0) {
1386  rv = oilrate / gasrate;
1387  }
1388  rv = std::min(rv, rvmax_perf[perf]);
1389 
1390  b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, 0.0 /*Rvw*/);
1391  }
1392  else {
1393  b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1394  }
1395  } else if (waterPresent) {
1396  //no oil
1397  const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers
1398  rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1399  if (waterrate > 0) {
1400  const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1401  double rvw = 0.0;
1402  if (gasrate > 0) {
1403  rvw = waterrate / gasrate;
1404  }
1405  rvw = std::min(rvw, rvwmax_perf[perf]);
1406 
1407  b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, 0.0 /*Rv*/, rvw);
1408  }
1409  else {
1410  b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1411  }
1412 
1413  } else {
1414  b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1415  }
1416  }
1417 
1418  if (oilPresent) {
1419  const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
1420  const int oilpos = oilCompIdx + perf * this->num_components_;
1421  if (gasPresent) {
1422  rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
1423  const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0);
1424  if (gasrate > 0) {
1425  const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
1426  double rs = 0.0;
1427  if (oilrate > 0) {
1428  rs = gasrate / oilrate;
1429  }
1430  rs = std::min(rs, rsmax_perf[perf]);
1431  b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
1432  } else {
1433  b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1434  }
1435  } else {
1436  b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1437  }
1438  }
1439 
1440  // Surface density.
1441  for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1442  if (!FluidSystem::phaseIsActive(phaseIdx)) {
1443  continue;
1444  }
1445 
1446  const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
1447  surf_dens_perf[this->num_components_ * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, fs.pvtRegionIndex() );
1448  }
1449 
1450  // We use cell values for solvent injector
1451  if constexpr (has_solvent) {
1452  b_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
1453  surf_dens_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventRefDensity();
1454  }
1455  }
1456  }
1457 
1458 
1459 
1460 
1461 
1462  template<typename TypeTag>
1463  ConvergenceReport
1465  getWellConvergence(const WellState& well_state,
1466  const std::vector<double>& B_avg,
1467  DeferredLogger& deferred_logger,
1468  const bool relax_tolerance) const
1469  {
1470  // the following implementation assume that the polymer is always after the w-o-g phases
1471  // For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells
1472  assert((int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
1473 
1474  std::vector<double> res;
1475  ConvergenceReport report = this->StdWellEval::getWellConvergence(well_state,
1476  B_avg,
1477  this->param_.max_residual_allowed_,
1478  this->param_.tolerance_wells_,
1479  this->param_.relaxed_tolerance_flow_well_,
1480  relax_tolerance,
1481  res,
1482  deferred_logger);
1483  checkConvergenceExtraEqs(res, report);
1484 
1485  return report;
1486  }
1487 
1488 
1489 
1490 
1491 
1492  template<typename TypeTag>
1493  void
1495  updateProductivityIndex(const Simulator& ebosSimulator,
1496  const WellProdIndexCalculator& wellPICalc,
1497  WellState& well_state,
1498  DeferredLogger& deferred_logger) const
1499  {
1500  auto fluidState = [&ebosSimulator, this](const int perf)
1501  {
1502  const auto cell_idx = this->well_cells_[perf];
1503  return ebosSimulator.model()
1504  .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
1505  };
1506 
1507  const int np = this->number_of_phases_;
1508  auto setToZero = [np](double* x) -> void
1509  {
1510  std::fill_n(x, np, 0.0);
1511  };
1512 
1513  auto addVector = [np](const double* src, double* dest) -> void
1514  {
1515  std::transform(src, src + np, dest, dest, std::plus<>{});
1516  };
1517 
1518  auto& ws = well_state.well(this->index_of_well_);
1519  auto& perf_data = ws.perf_data;
1520  auto* wellPI = ws.productivity_index.data();
1521  auto* connPI = perf_data.prod_index.data();
1522 
1523  setToZero(wellPI);
1524 
1525  const auto preferred_phase = this->well_ecl_.getPreferredPhase();
1526  auto subsetPerfID = 0;
1527 
1528  for (const auto& perf : *this->perf_data_) {
1529  auto allPerfID = perf.ecl_index;
1530 
1531  auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
1532  {
1533  return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
1534  };
1535 
1536  std::vector<EvalWell> mob(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.0});
1537  getMobilityEval(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
1538 
1539  const auto& fs = fluidState(subsetPerfID);
1540  setToZero(connPI);
1541 
1542  if (this->isInjector()) {
1543  this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
1544  mob, connPI, deferred_logger);
1545  }
1546  else { // Production or zero flow rate
1547  this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
1548  }
1549 
1550  addVector(connPI, wellPI);
1551 
1552  ++subsetPerfID;
1553  connPI += np;
1554  }
1555 
1556  // Sum with communication in case of distributed well.
1557  const auto& comm = this->parallel_well_info_.communication();
1558  if (comm.size() > 1) {
1559  comm.sum(wellPI, np);
1560  }
1561 
1562  assert ((static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
1563  "Internal logic error in processing connections for PI/II");
1564  }
1565 
1566 
1567 
1568  template<typename TypeTag>
1569  void
1570  StandardWell<TypeTag>::
1571  computeWellConnectionDensitesPressures(const Simulator& ebosSimulator,
1572  const WellState& well_state,
1573  const std::vector<double>& b_perf,
1574  const std::vector<double>& rsmax_perf,
1575  const std::vector<double>& rvmax_perf,
1576  const std::vector<double>& rvwmax_perf,
1577  const std::vector<double>& surf_dens_perf,
1578  DeferredLogger& deferred_logger)
1579  {
1580  // Compute densities
1581  const int nperf = this->number_of_perforations_;
1582  const int np = this->number_of_phases_;
1583  std::vector<double> perfRates(b_perf.size(),0.0);
1584  const auto& ws = well_state.well(this->index_of_well_);
1585  const auto& perf_data = ws.perf_data;
1586  const auto& perf_rates_state = perf_data.phase_rates;
1587 
1588  for (int perf = 0; perf < nperf; ++perf) {
1589  for (int comp = 0; comp < np; ++comp) {
1590  perfRates[perf * this->num_components_ + comp] = perf_rates_state[perf * np + this->ebosCompIdxToFlowCompIdx(comp)];
1591  }
1592  }
1593 
1594  if constexpr (has_solvent) {
1595  const auto& solvent_perf_rates_state = perf_data.solvent_rates;
1596  for (int perf = 0; perf < nperf; ++perf) {
1597  perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = solvent_perf_rates_state[perf];
1598  }
1599  }
1600 
1601  // for producers where all perforations have zero rate we
1602  // approximate the perforation mixture using the mobility ratio
1603  // and weight the perforations using the well transmissibility.
1604  bool all_zero = std::all_of(perfRates.begin(), perfRates.end(), [](double val) { return val == 0.0; });
1605  const auto& comm = this->parallel_well_info_.communication();
1606  if (comm.size() > 1)
1607  {
1608  all_zero = (comm.min(all_zero ? 1 : 0) == 1);
1609  }
1610 
1611  if ( all_zero && this->isProducer() ) {
1612  double total_tw = 0;
1613  for (int perf = 0; perf < nperf; ++perf) {
1614  total_tw += this->well_index_[perf];
1615  }
1616  if (comm.size() > 1)
1617  {
1618  total_tw = comm.sum(total_tw);
1619  }
1620  for (int perf = 0; perf < nperf; ++perf) {
1621  const int cell_idx = this->well_cells_[perf];
1622  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1623  const auto& fs = intQuants.fluidState();
1624  const double well_tw_fraction = this->well_index_[perf] / total_tw;
1625  double total_mobility = 0.0;
1626  for (int p = 0; p < np; ++p) {
1627  int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1628  total_mobility += fs.invB(ebosPhaseIdx).value() * intQuants.mobility(ebosPhaseIdx).value();
1629  }
1630  if constexpr (has_solvent) {
1631  total_mobility += intQuants.solventInverseFormationVolumeFactor().value() * intQuants.solventMobility().value();
1632  }
1633  for (int p = 0; p < np; ++p) {
1634  int ebosPhaseIdx = this->flowPhaseToEbosPhaseIdx(p);
1635  perfRates[perf * this->num_components_ + p] = well_tw_fraction * intQuants.mobility(ebosPhaseIdx).value() / total_mobility;
1636  }
1637  if constexpr (has_solvent) {
1638  perfRates[perf * this->num_components_ + Indices::contiSolventEqIdx] = well_tw_fraction * intQuants.solventInverseFormationVolumeFactor().value() / total_mobility;
1639  }
1640  }
1641  }
1642 
1643  this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
1644 
1645  this->computeConnectionPressureDelta();
1646  }
1647 
1648 
1649 
1650 
1651 
1652  template<typename TypeTag>
1653  void
1654  StandardWell<TypeTag>::
1655  computeWellConnectionPressures(const Simulator& ebosSimulator,
1656  const WellState& well_state,
1657  DeferredLogger& deferred_logger)
1658  {
1659  // 1. Compute properties required by computeConnectionPressureDelta().
1660  // Note that some of the complexity of this part is due to the function
1661  // taking std::vector<double> arguments, and not Eigen objects.
1662  std::vector<double> b_perf;
1663  std::vector<double> rsmax_perf;
1664  std::vector<double> rvmax_perf;
1665  std::vector<double> rvwmax_perf;
1666  std::vector<double> surf_dens_perf;
1667  computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf);
1668  computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, rvwmax_perf, surf_dens_perf, deferred_logger);
1669  }
1670 
1671 
1672 
1673 
1674 
1675  template<typename TypeTag>
1676  void
1677  StandardWell<TypeTag>::
1678  solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
1679  {
1680  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1681 
1682  // We assemble the well equations, then we check the convergence,
1683  // which is why we do not put the assembleWellEq here.
1684  BVectorWell dx_well(1);
1685  dx_well[0].resize(this->numWellEq_);
1686  this->invDuneD_.mv(this->resWell_, dx_well);
1687 
1688  updateWellState(dx_well, well_state, deferred_logger);
1689  }
1690 
1691 
1692 
1693 
1694 
1695  template<typename TypeTag>
1696  void
1697  StandardWell<TypeTag>::
1698  calculateExplicitQuantities(const Simulator& ebosSimulator,
1699  const WellState& well_state,
1700  DeferredLogger& deferred_logger)
1701  {
1702  updatePrimaryVariables(well_state, deferred_logger);
1703  initPrimaryVariablesEvaluation();
1704  computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
1705  this->computeAccumWell();
1706  }
1707 
1708 
1709 
1710  template<typename TypeTag>
1711  void
1713  apply(const BVector& x, BVector& Ax) const
1714  {
1715  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1716 
1717  if (this->param_.matrix_add_well_contributions_)
1718  {
1719  // Contributions are already in the matrix itself
1720  return;
1721  }
1722  assert( this->Bx_.size() == this->duneB_.N() );
1723  assert( this->invDrw_.size() == this->invDuneD_.N() );
1724 
1725  // Bx_ = duneB_ * x
1726  this->parallelB_.mv(x, this->Bx_);
1727 
1728  // invDBx = invDuneD_ * Bx_
1729  // TODO: with this, we modified the content of the invDrw_.
1730  // Is it necessary to do this to save some memory?
1731  BVectorWell& invDBx = this->invDrw_;
1732  this->invDuneD_.mv(this->Bx_, invDBx);
1733 
1734  // Ax = Ax - duneC_^T * invDBx
1735  this->duneC_.mmtv(invDBx,Ax);
1736  }
1737 
1738 
1739 
1740 
1741  template<typename TypeTag>
1742  void
1744  apply(BVector& r) const
1745  {
1746  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1747 
1748  assert( this->invDrw_.size() == this->invDuneD_.N() );
1749 
1750  // invDrw_ = invDuneD_ * resWell_
1751  this->invDuneD_.mv(this->resWell_, this->invDrw_);
1752  // r = r - duneC_^T * invDrw_
1753  this->duneC_.mmtv(this->invDrw_, r);
1754  }
1755 
1756  template<typename TypeTag>
1757  void
1759  recoverSolutionWell(const BVector& x, BVectorWell& xw) const
1760  {
1761  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1762 
1763  BVectorWell resWell = this->resWell_;
1764  // resWell = resWell - B * x
1765  this->parallelB_.mmv(x, resWell);
1766  // xw = D^-1 * resWell
1767  this->invDuneD_.mv(resWell, xw);
1768  }
1769 
1770 
1771 
1772 
1773 
1774  template<typename TypeTag>
1775  void
1778  WellState& well_state,
1779  DeferredLogger& deferred_logger) const
1780  {
1781  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
1782 
1783  BVectorWell xw(1);
1784  xw[0].resize(this->numWellEq_);
1785 
1786  recoverSolutionWell(x, xw);
1787  updateWellState(xw, well_state, deferred_logger);
1788  }
1789 
1790 
1791 
1792 
1793  template<typename TypeTag>
1794  void
1796  computeWellRatesWithBhp(const Simulator& ebosSimulator,
1797  const double& bhp,
1798  std::vector<double>& well_flux,
1799  DeferredLogger& deferred_logger) const
1800  {
1801 
1802  const int np = this->number_of_phases_;
1803  well_flux.resize(np, 0.0);
1804 
1805  const bool allow_cf = this->getAllowCrossFlow();
1806 
1807  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
1808  const int cell_idx = this->well_cells_[perf];
1809  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1810  // flux for each perforation
1811  std::vector<Scalar> mob(this->num_components_, 0.);
1812  getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
1813  double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
1814  const double Tw = this->well_index_[perf] * trans_mult;
1815 
1816  std::vector<Scalar> cq_s(this->num_components_, 0.);
1817  computePerfRateScalar(intQuants, mob, bhp, Tw, perf, allow_cf,
1818  cq_s, deferred_logger);
1819 
1820  for(int p = 0; p < np; ++p) {
1821  well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p];
1822  }
1823  }
1824  this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size());
1825  }
1826 
1827 
1828 
1829  template<typename TypeTag>
1830  void
1831  StandardWell<TypeTag>::
1832  computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
1833  const double& bhp,
1834  std::vector<double>& well_flux,
1835  DeferredLogger& deferred_logger) const
1836  {
1837 
1838  // iterate to get a more accurate well density
1839  // create a copy of the well_state to use. If the operability checking is sucessful, we use this one
1840  // to replace the original one
1841  WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
1842  const auto& group_state = ebosSimulator.problem().wellModel().groupState();
1843  auto& ws = well_state_copy.well(this->index_of_well_);
1844 
1845  // Set current control to bhp, and bhp value in state, modify bhp limit in control object.
1846  if (this->well_ecl_.isInjector()) {
1847  ws.injection_cmode = Well::InjectorCMode::BHP;
1848  } else {
1849  ws.production_cmode = Well::ProducerCMode::BHP;
1850  }
1851  ws.bhp = bhp;
1852 
1853  // initialized the well rates with the potentials i.e. the well rates based on bhp
1854  const int np = this->number_of_phases_;
1855  const double sign = this->well_ecl_.isInjector() ? 1.0 : -1.0;
1856  for (int phase = 0; phase < np; ++phase){
1857  well_state_copy.wellRates(this->index_of_well_)[phase]
1858  = sign * ws.well_potentials[phase];
1859  }
1860  // creating a copy of the well itself, to avoid messing up the explicit informations
1861  // during this copy, the only information not copied properly is the well controls
1862  StandardWell<TypeTag> well(*this);
1863  well.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
1864 
1865  const double dt = ebosSimulator.timeStepSize();
1866  bool converged = well.iterateWellEquations(ebosSimulator, dt, well_state_copy, group_state, deferred_logger);
1867  if (!converged) {
1868  const std::string msg = " well " + name() + " did not get converged during well potential calculations "
1869  " potentials are computed based on unconverged solution";
1870  deferred_logger.debug(msg);
1871  }
1872  well.updatePrimaryVariables(well_state_copy, deferred_logger);
1873  well.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
1874  well.initPrimaryVariablesEvaluation();
1875  well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
1876  }
1877 
1878 
1879 
1880 
1881  template<typename TypeTag>
1882  std::vector<double>
1883  StandardWell<TypeTag>::
1884  computeWellPotentialWithTHP(const Simulator& ebos_simulator,
1885  DeferredLogger& deferred_logger,
1886  const WellState &well_state) const
1887  {
1888  std::vector<double> potentials(this->number_of_phases_, 0.0);
1889  const auto& summary_state = ebos_simulator.vanguard().summaryState();
1890 
1891  const auto& well = this->well_ecl_;
1892  if (well.isInjector()){
1893  const auto& controls = this->well_ecl_.injectionControls(summary_state);
1894  auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
1895  if (bhp_at_thp_limit) {
1896  const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
1897  computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1898  } else {
1899  deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
1900  "Failed in getting converged thp based potential calculation for well "
1901  + name() + ". Instead the bhp based value is used");
1902  const double bhp = controls.bhp_limit;
1903  computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1904  }
1905  } else {
1906  computeWellRatesWithThpAlqProd(
1907  ebos_simulator, summary_state,
1908  deferred_logger, potentials, this->getALQ(well_state)
1909  );
1910  }
1911 
1912  return potentials;
1913  }
1914 
1915  template<typename TypeTag>
1916  double
1917  StandardWell<TypeTag>::
1918  computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
1919  const SummaryState &summary_state,
1920  DeferredLogger &deferred_logger,
1921  std::vector<double> &potentials,
1922  double alq) const
1923  {
1924  double bhp;
1925  auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
1926  ebos_simulator, summary_state, deferred_logger, alq);
1927  if (bhp_at_thp_limit) {
1928  const auto& controls = this->well_ecl_.productionControls(summary_state);
1929  bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
1930  computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1931  }
1932  else {
1933  deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
1934  "Failed in getting converged thp based potential calculation for well "
1935  + name() + ". Instead the bhp based value is used");
1936  const auto& controls = this->well_ecl_.productionControls(summary_state);
1937  bhp = controls.bhp_limit;
1938  computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
1939  }
1940  return bhp;
1941  }
1942 
1943  template<typename TypeTag>
1944  void
1945  StandardWell<TypeTag>::
1946  computeWellRatesWithThpAlqProd(const Simulator &ebos_simulator,
1947  const SummaryState &summary_state,
1948  DeferredLogger &deferred_logger,
1949  std::vector<double> &potentials,
1950  double alq) const
1951  {
1952  /*double bhp =*/
1953  computeWellRatesAndBhpWithThpAlqProd(ebos_simulator,
1954  summary_state,
1955  deferred_logger,
1956  potentials,
1957  alq);
1958  }
1959 
1960  template<typename TypeTag>
1961  void
1963  computeWellPotentials(const Simulator& ebosSimulator,
1964  const WellState& well_state,
1965  std::vector<double>& well_potentials,
1966  DeferredLogger& deferred_logger) // const
1967  {
1968  const int np = this->number_of_phases_;
1969  well_potentials.resize(np, 0.0);
1970 
1971  if (this->wellIsStopped()) {
1972  return;
1973  }
1974 
1975  this->operability_status_.has_negative_potentials = false;
1976  // If the well is pressure controlled the potential equals the rate.
1977  bool thp_controlled_well = false;
1978  bool bhp_controlled_well = false;
1979  const auto& ws = well_state.well(this->index_of_well_);
1980  if (this->isInjector()) {
1981  const Well::InjectorCMode& current = ws.injection_cmode;
1982  if (current == Well::InjectorCMode::THP) {
1983  thp_controlled_well = true;
1984  }
1985  if (current == Well::InjectorCMode::BHP) {
1986  bhp_controlled_well = true;
1987  }
1988  } else {
1989  const Well::ProducerCMode& current = ws.production_cmode;
1990  if (current == Well::ProducerCMode::THP) {
1991  thp_controlled_well = true;
1992  }
1993  if (current == Well::ProducerCMode::BHP) {
1994  bhp_controlled_well = true;
1995  }
1996  }
1997  if (!this->changed_to_open_this_step_ && (thp_controlled_well || bhp_controlled_well)) {
1998 
1999  double total_rate = 0.0;
2000  const double sign = this->isInjector() ? 1.0:-1.0;
2001  for (int phase = 0; phase < np; ++phase){
2002  total_rate += sign * ws.surface_rates[phase];
2003  }
2004  // for pressure controlled wells the well rates are the potentials
2005  // if the rates are trivial we are most probably looking at the newly
2006  // opened well and we therefore make the affort of computing the potentials anyway.
2007  if (total_rate > 0) {
2008  for (int phase = 0; phase < np; ++phase){
2009  well_potentials[phase] = sign * ws.surface_rates[phase];
2010  }
2011  return;
2012  }
2013  }
2014 
2015  // does the well have a THP related constraint?
2016  const auto& summaryState = ebosSimulator.vanguard().summaryState();
2017  if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
2018  // get the bhp value based on the bhp constraints
2019  double bhp = this->mostStrictBhpFromBhpLimits(summaryState);
2020 
2021  // In some very special cases the bhp pressure target are
2022  // temporary violated. This may lead to too small or negative potentials
2023  // that could lead to premature shutting of wells.
2024  // As a remedy the bhp that gives the largest potential is used.
2025  // For converged cases, ws.bhp <=bhp for injectors and ws.bhp >= bhp,
2026  // and the potentials will be computed using the limit as expected.
2027  if (this->isInjector())
2028  bhp = std::max(ws.bhp, bhp);
2029  else
2030  bhp = std::min(ws.bhp, bhp);
2031 
2032  assert(std::abs(bhp) != std::numeric_limits<double>::max());
2033  computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger);
2034  } else {
2035  // the well has a THP related constraint
2036  well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
2037  }
2038 
2039  const double sign = this->isInjector() ? 1.0:-1.0;
2040  double total_potential = 0.0;
2041  for (int phase = 0; phase < np; ++phase){
2042  well_potentials[phase] *= sign;
2043  total_potential += well_potentials[phase];
2044  }
2045  if (total_potential < 0.0 && this->param_.check_well_operability_) {
2046  // wells with negative potentials are not operable
2047  this->operability_status_.has_negative_potentials = true;
2048  const std::string msg = std::string("well ") + this->name() + std::string(": has negative potentials and is not operable");
2049  deferred_logger.warning("NEGATIVE_POTENTIALS_INOPERABLE", msg);
2050  }
2051  }
2052 
2053 
2054 
2055 
2056 
2057  template<typename TypeTag>
2058  void
2060  updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const
2061  {
2062  this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger);
2063  if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
2064 
2065  // other primary variables related to polymer injection
2066  if constexpr (Base::has_polymermw) {
2067  if (this->isInjector()) {
2068  const auto& ws = well_state.well(this->index_of_well_);
2069  const auto& perf_data = ws.perf_data;
2070  const auto& water_velocity = perf_data.water_velocity;
2071  const auto& skin_pressure = perf_data.skin_pressure;
2072  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
2073  this->primary_variables_[Bhp + 1 + perf] = water_velocity[perf];
2074  this->primary_variables_[Bhp + 1 + this->number_of_perforations_ + perf] = skin_pressure[perf];
2075  }
2076  }
2077  }
2078  for (double v : this->primary_variables_) {
2079  if(!isfinite(v))
2080  OPM_DEFLOG_THROW(NumericalIssue, "Infinite primary variable after update from wellState well: " << this->name(), deferred_logger);
2081  }
2082  }
2083 
2084 
2085 
2086 
2087  template<typename TypeTag>
2088  double
2089  StandardWell<TypeTag>::
2090  getRefDensity() const
2091  {
2092  return this->perf_densities_[0];
2093  }
2094 
2095 
2096 
2097 
2098  template<typename TypeTag>
2099  void
2100  StandardWell<TypeTag>::
2101  updateWaterMobilityWithPolymer(const Simulator& ebos_simulator,
2102  const int perf,
2103  std::vector<EvalWell>& mob,
2104  DeferredLogger& deferred_logger) const
2105  {
2106  const int cell_idx = this->well_cells_[perf];
2107  const auto& int_quant = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
2108  const EvalWell polymer_concentration = this->extendEval(int_quant.polymerConcentration());
2109 
2110  // TODO: not sure should based on the well type or injecting/producing peforations
2111  // it can be different for crossflow
2112  if (this->isInjector()) {
2113  // assume fully mixing within injecting wellbore
2114  const auto& visc_mult_table = PolymerModule::plyviscViscosityMultiplierTable(int_quant.pvtRegionIndex());
2115  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2116  mob[waterCompIdx] /= (this->extendEval(int_quant.waterViscosityCorrection()) * visc_mult_table.eval(polymer_concentration, /*extrapolate=*/true) );
2117  }
2118 
2119  if (PolymerModule::hasPlyshlog()) {
2120  // we do not calculate the shear effects for injection wells when they do not
2121  // inject polymer.
2122  if (this->isInjector() && this->wpolymer() == 0.) {
2123  return;
2124  }
2125  // compute the well water velocity with out shear effects.
2126  // TODO: do we need to turn on crossflow here?
2127  const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebos_simulator);
2128  const EvalWell& bhp = this->getBhp();
2129 
2130  std::vector<EvalWell> cq_s(this->num_components_, {this->numWellEq_ + Indices::numEq, 0.});
2131  double perf_dis_gas_rate = 0.;
2132  double perf_vap_oil_rate = 0.;
2133  double perf_vap_wat_rate = 0.;
2134  double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quant, cell_idx);
2135  const double Tw = this->well_index_[perf] * trans_mult;
2136  computePerfRateEval(int_quant, mob, bhp, Tw, perf, allow_cf,
2137  cq_s, perf_dis_gas_rate, perf_vap_oil_rate, perf_vap_wat_rate, deferred_logger);
2138  // TODO: make area a member
2139  const double area = 2 * M_PI * this->perf_rep_radius_[perf] * this->perf_length_[perf];
2140  const auto& material_law_manager = ebos_simulator.problem().materialLawManager();
2141  const auto& scaled_drainage_info =
2142  material_law_manager->oilWaterScaledEpsInfoDrainage(cell_idx);
2143  const double swcr = scaled_drainage_info.Swcr;
2144  const EvalWell poro = this->extendEval(int_quant.porosity());
2145  const EvalWell sw = this->extendEval(int_quant.fluidState().saturation(FluidSystem::waterPhaseIdx));
2146  // guard against zero porosity and no water
2147  const EvalWell denom = max( (area * poro * (sw - swcr)), 1e-12);
2148  const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2149  EvalWell water_velocity = cq_s[waterCompIdx] / denom * this->extendEval(int_quant.fluidState().invB(FluidSystem::waterPhaseIdx));
2150 
2151  if (PolymerModule::hasShrate()) {
2152  // the equation for the water velocity conversion for the wells and reservoir are from different version
2153  // of implementation. It can be changed to be more consistent when possible.
2154  water_velocity *= PolymerModule::shrate( int_quant.pvtRegionIndex() ) / this->bore_diameters_[perf];
2155  }
2156  const EvalWell shear_factor = PolymerModule::computeShearFactor(polymer_concentration,
2157  int_quant.pvtRegionIndex(),
2158  water_velocity);
2159  // modify the mobility with the shear factor.
2160  mob[waterCompIdx] /= shear_factor;
2161  }
2162  }
2163 
2164  template<typename TypeTag>
2165  void
2166  StandardWell<TypeTag>::addWellContributions(SparseMatrixAdapter& jacobian) const
2167  {
2168  // We need to change matrx A as follows
2169  // A -= C^T D^-1 B
2170  // D is diagonal
2171  // B and C have 1 row, nc colums and nonzero
2172  // at (0,j) only if this well has a perforation at cell j.
2173  typename SparseMatrixAdapter::MatrixBlock tmpMat;
2174  Dune::DynamicMatrix<Scalar> tmp;
2175  for ( auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC )
2176  {
2177  const auto row_index = colC.index();
2178 
2179  for ( auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB )
2180  {
2181  detail::multMatrix(this->invDuneD_[0][0], (*colB), tmp);
2182  detail::negativeMultMatrixTransposed((*colC), tmp, tmpMat);
2183  jacobian.addToBlock( row_index, colB.index(), tmpMat );
2184  }
2185  }
2186  }
2187 
2188 
2189  template <typename TypeTag>
2190  void
2191  StandardWell<TypeTag>::addWellPressureEquations(PressureMatrix& jacobian,
2192  const BVector& weights,
2193  const int pressureVarIndex,
2194  const bool use_well_weights,
2195  const WellState& well_state) const
2196  {
2197  // This adds pressure quation for cpr
2198  // For use_well_weights=true
2199  // weights lamda = inv(D)'v v = 0 v(bhpInd) = 1
2200  // the well equations are summed i lambda' B(:,pressureVarINd) -> B lambda'*D(:,bhpInd) -> D
2201  // For use_well_weights = false
2202  // weights lambda = \sum_i w /n where ths sum is over weights of all perforation cells
2203  // in the case of pressure controlled trivial equations are used and bhp C=B=0
2204  // then the flow part of the well equations are summed lambda'*B(1:n,pressureVarInd) -> B lambda'*D(1:n,bhpInd) -> D
2205  // For bouth
2206  // C -> w'C(:,bhpInd) where w is weights of the perforation cell
2207 
2208  // Add the well contributions in cpr
2209  // use_well_weights is a quasiimpes formulation which is not implemented in multisegment
2210  int bhp_var_index = Bhp;
2211  int nperf = 0;
2212  auto cell_weights = weights[0];// not need for not(use_well_weights)
2213  cell_weights = 0.0;
2214  assert(this->duneC_.M() == weights.size());
2215  const int welldof_ind = this->duneC_.M() + this->index_of_well_;
2216  // do not assume anything about pressure controlled with use_well_weights (work fine with the assumtion also)
2217  if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
2218  // make coupling for reservoir to well
2219  for (auto colC = this->duneC_[0].begin(), endC = this->duneC_[0].end(); colC != endC; ++colC) {
2220  const auto row_ind = colC.index();
2221  const auto& bw = weights[row_ind];
2222  double matel = 0;
2223  assert((*colC).M() == bw.size());
2224  for (size_t i = 0; i < bw.size(); ++i) {
2225  matel += (*colC)[bhp_var_index][i] * bw[i];
2226  }
2227 
2228  jacobian[row_ind][welldof_ind] = matel;
2229  cell_weights += bw;
2230  nperf += 1;
2231  }
2232  }
2233  cell_weights /= nperf;
2234 
2235  BVectorWell bweights(1);
2236  size_t blockSz = this->numWellEq_;
2237  bweights[0].resize(blockSz);
2238  bweights[0] = 0.0;
2239  double diagElem = 0;
2240  {
2241  if ( use_well_weights ){
2242  // calculate weighs and set diagonal element
2243  //NB! use this options without treating pressure controlled separated
2244  //NB! calculate quasiimpes well weights NB do not work well with trueimpes reservoir weights
2245  double abs_max = 0;
2246  BVectorWell rhs(1);
2247  rhs[0].resize(blockSz);
2248  rhs[0][bhp_var_index] = 1.0;
2249  DiagMatrixBlockWellType inv_diag_block = this->invDuneD_[0][0];
2250  DiagMatrixBlockWellType inv_diag_block_transpose = Opm::wellhelpers::transposeDenseDynMatrix(inv_diag_block);
2251  for (size_t i = 0; i < blockSz; ++i) {
2252  bweights[0][i] = 0;
2253  for (size_t j = 0; j < blockSz; ++j) {
2254  bweights[0][i] += inv_diag_block_transpose[i][j]*rhs[0][j];
2255  }
2256  abs_max = std::max(abs_max, std::fabs(bweights[0][i]));
2257  }
2258  assert( abs_max > 0.0 );
2259  for (size_t i = 0; i < blockSz; ++i) {
2260  bweights[0][i] /= abs_max;
2261  }
2262  diagElem = 1.0/abs_max;
2263  }else{
2264  // set diagonal element
2265  if( this->isPressureControlled(well_state) ){
2266  bweights[0][blockSz-1] = 1.0;
2267  diagElem = 1.0;// better scaling could have used the calculation below if weights were calculated
2268  }else{
2269  for (size_t i = 0; i < cell_weights.size(); ++i) {
2270  bweights[0][i] = cell_weights[i];
2271  }
2272  bweights[0][blockSz-1] = 0.0;
2273  diagElem = 0.0;
2274  const auto& locmat = this->duneD_[0][0];
2275  for (size_t i = 0; i < cell_weights.size(); ++i) {
2276  diagElem += locmat[i][bhp_var_index]*cell_weights[i];
2277  }
2278 
2279  }
2280  }
2281  }
2282  //
2283  jacobian[welldof_ind][welldof_ind] = diagElem;
2284  // set the matrix elements for well reservoir coupling
2285  if( not( this->isPressureControlled(well_state) ) || use_well_weights ){
2286  for (auto colB = this->duneB_[0].begin(), endB = this->duneB_[0].end(); colB != endB; ++colB) {
2287  const auto col_index = colB.index();
2288  const auto& bw = bweights[0];
2289  double matel = 0;
2290  for (size_t i = 0; i < bw.size(); ++i) {
2291  matel += (*colB)[i][pressureVarIndex] * bw[i];
2292  }
2293  jacobian[welldof_ind][col_index] = matel;
2294  }
2295  }
2296  }
2297 
2298 
2299 
2300  template<typename TypeTag>
2301  typename StandardWell<TypeTag>::EvalWell
2302  StandardWell<TypeTag>::
2303  pskinwater(const double throughput,
2304  const EvalWell& water_velocity,
2305  DeferredLogger& deferred_logger) const
2306  {
2307  if constexpr (Base::has_polymermw) {
2308  const int water_table_id = this->well_ecl_.getPolymerProperties().m_skprwattable;
2309  if (water_table_id <= 0) {
2310  OPM_DEFLOG_THROW(std::runtime_error, "Unused SKPRWAT table id used for well " << name(), deferred_logger);
2311  }
2312  const auto& water_table_func = PolymerModule::getSkprwatTable(water_table_id);
2313  const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2314  // the skin pressure when injecting water, which also means the polymer concentration is zero
2315  EvalWell pskin_water(this->numWellEq_ + Indices::numEq, 0.0);
2316  pskin_water = water_table_func.eval(throughput_eval, water_velocity);
2317  return pskin_water;
2318  } else {
2319  OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
2320  "while injecting skin pressure is requested for well " << name(), deferred_logger);
2321  }
2322  }
2323 
2324 
2325 
2326 
2327 
2328  template<typename TypeTag>
2329  typename StandardWell<TypeTag>::EvalWell
2330  StandardWell<TypeTag>::
2331  pskin(const double throughput,
2332  const EvalWell& water_velocity,
2333  const EvalWell& poly_inj_conc,
2334  DeferredLogger& deferred_logger) const
2335  {
2336  if constexpr (Base::has_polymermw) {
2337  const double sign = water_velocity >= 0. ? 1.0 : -1.0;
2338  const EvalWell water_velocity_abs = abs(water_velocity);
2339  if (poly_inj_conc == 0.) {
2340  return sign * pskinwater(throughput, water_velocity_abs, deferred_logger);
2341  }
2342  const int polymer_table_id = this->well_ecl_.getPolymerProperties().m_skprpolytable;
2343  if (polymer_table_id <= 0) {
2344  OPM_DEFLOG_THROW(std::runtime_error, "Unavailable SKPRPOLY table id used for well " << name(), deferred_logger);
2345  }
2346  const auto& skprpolytable = PolymerModule::getSkprpolyTable(polymer_table_id);
2347  const double reference_concentration = skprpolytable.refConcentration;
2348  const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2349  // the skin pressure when injecting water, which also means the polymer concentration is zero
2350  EvalWell pskin_poly(this->numWellEq_ + Indices::numEq, 0.0);
2351  pskin_poly = skprpolytable.table_func.eval(throughput_eval, water_velocity_abs);
2352  if (poly_inj_conc == reference_concentration) {
2353  return sign * pskin_poly;
2354  }
2355  // poly_inj_conc != reference concentration of the table, then some interpolation will be required
2356  const EvalWell pskin_water = pskinwater(throughput, water_velocity_abs, deferred_logger);
2357  const EvalWell pskin = pskin_water + (pskin_poly - pskin_water) / reference_concentration * poly_inj_conc;
2358  return sign * pskin;
2359  } else {
2360  OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
2361  "while injecting skin pressure is requested for well " << name(), deferred_logger);
2362  }
2363  }
2364 
2365 
2366 
2367 
2368 
2369  template<typename TypeTag>
2370  typename StandardWell<TypeTag>::EvalWell
2371  StandardWell<TypeTag>::
2372  wpolymermw(const double throughput,
2373  const EvalWell& water_velocity,
2374  DeferredLogger& deferred_logger) const
2375  {
2376  if constexpr (Base::has_polymermw) {
2377  const int table_id = this->well_ecl_.getPolymerProperties().m_plymwinjtable;
2378  const auto& table_func = PolymerModule::getPlymwinjTable(table_id);
2379  const EvalWell throughput_eval(this->numWellEq_ + Indices::numEq, throughput);
2380  EvalWell molecular_weight(this->numWellEq_ + Indices::numEq, 0.);
2381  if (this->wpolymer() == 0.) { // not injecting polymer
2382  return molecular_weight;
2383  }
2384  molecular_weight = table_func.eval(throughput_eval, abs(water_velocity));
2385  return molecular_weight;
2386  } else {
2387  OPM_DEFLOG_THROW(std::runtime_error, "Polymermw is not activated, "
2388  "while injecting polymer molecular weight is requested for well " << name(), deferred_logger);
2389  }
2390  }
2391 
2392 
2393 
2394 
2395 
2396  template<typename TypeTag>
2397  void
2398  StandardWell<TypeTag>::
2399  updateWaterThroughput(const double dt, WellState &well_state) const
2400  {
2401  if constexpr (Base::has_polymermw) {
2402  if (this->isInjector()) {
2403  auto& ws = well_state.well(this->index_of_well_);
2404  auto& perf_water_throughput = ws.perf_data.water_throughput;
2405  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
2406  const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
2407  // we do not consider the formation damage due to water flowing from reservoir into wellbore
2408  if (perf_water_vel > 0.) {
2409  perf_water_throughput[perf] += perf_water_vel * dt;
2410  }
2411  }
2412  }
2413  }
2414  }
2415 
2416 
2417 
2418 
2419 
2420  template<typename TypeTag>
2421  void
2422  StandardWell<TypeTag>::
2423  handleInjectivityRate(const Simulator& ebosSimulator,
2424  const int perf,
2425  std::vector<EvalWell>& cq_s) const
2426  {
2427  const int cell_idx = this->well_cells_[perf];
2428  const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
2429  const auto& fs = int_quants.fluidState();
2430  const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2431  const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2432  const int wat_vel_index = Bhp + 1 + perf;
2433  const unsigned water_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
2434 
2435  // water rate is update to use the form from water velocity, since water velocity is
2436  // a primary variable now
2437  cq_s[water_comp_idx] = area * this->primary_variables_evaluation_[wat_vel_index] * b_w;
2438  }
2439 
2440 
2441 
2442 
2443  template<typename TypeTag>
2444  void
2445  StandardWell<TypeTag>::
2446  handleInjectivityEquations(const Simulator& ebosSimulator,
2447  const WellState& well_state,
2448  const int perf,
2449  const EvalWell& water_flux_s,
2450  DeferredLogger& deferred_logger)
2451  {
2452  const int cell_idx = this->well_cells_[perf];
2453  const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
2454  const auto& fs = int_quants.fluidState();
2455  const EvalWell b_w = this->extendEval(fs.invB(FluidSystem::waterPhaseIdx));
2456  const EvalWell water_flux_r = water_flux_s / b_w;
2457  const double area = M_PI * this->bore_diameters_[perf] * this->perf_length_[perf];
2458  const EvalWell water_velocity = water_flux_r / area;
2459  const int wat_vel_index = Bhp + 1 + perf;
2460 
2461  // equation for the water velocity
2462  const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
2463  this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
2464 
2465  const auto& ws = well_state.well(this->index_of_well_);
2466  const auto& perf_data = ws.perf_data;
2467  const auto& perf_water_throughput = perf_data.water_throughput;
2468  const double throughput = perf_water_throughput[perf];
2469  const int pskin_index = Bhp + 1 + this->number_of_perforations_ + perf;
2470 
2471  EvalWell poly_conc(this->numWellEq_ + Indices::numEq, 0.0);
2472  poly_conc.setValue(this->wpolymer());
2473 
2474  // equation for the skin pressure
2475  const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
2476  - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
2477 
2478  this->resWell_[0][pskin_index] = eq_pskin.value();
2479  for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
2480  this->duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
2481  this->duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
2482  }
2483 
2484  // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
2485  for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
2486  this->duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
2487  }
2488  }
2489 
2490 
2491 
2492 
2493 
2494  template<typename TypeTag>
2495  void
2496  StandardWell<TypeTag>::
2497  checkConvergenceExtraEqs(const std::vector<double>& res,
2498  ConvergenceReport& report) const
2499  {
2500  // if different types of extra equations are involved, this function needs to be refactored further
2501 
2502  // checking the convergence of the extra equations related to polymer injectivity
2503  if constexpr (Base::has_polymermw) {
2504  this->checkConvergencePolyMW(res, report, this->param_.max_residual_allowed_);
2505  }
2506  }
2507 
2508 
2509 
2510 
2511 
2512  template<typename TypeTag>
2513  void
2514  StandardWell<TypeTag>::
2515  updateConnectionRatePolyMW(const EvalWell& cq_s_poly,
2516  const IntensiveQuantities& int_quants,
2517  const WellState& well_state,
2518  const int perf,
2519  std::vector<RateVector>& connectionRates,
2520  DeferredLogger& deferred_logger) const
2521  {
2522  // the source term related to transport of molecular weight
2523  EvalWell cq_s_polymw = cq_s_poly;
2524  if (this->isInjector()) {
2525  const int wat_vel_index = Bhp + 1 + perf;
2526  const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
2527  if (water_velocity > 0.) { // injecting
2528  const auto& ws = well_state.well(this->index_of_well_);
2529  const auto& perf_water_throughput = ws.perf_data.water_throughput;
2530  const double throughput = perf_water_throughput[perf];
2531  const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
2532  cq_s_polymw *= molecular_weight;
2533  } else {
2534  // we do not consider the molecular weight from the polymer
2535  // going-back to the wellbore through injector
2536  cq_s_polymw *= 0.;
2537  }
2538  } else if (this->isProducer()) {
2539  if (cq_s_polymw < 0.) {
2540  cq_s_polymw *= this->extendEval(int_quants.polymerMoleWeight() );
2541  } else {
2542  // we do not consider the molecular weight from the polymer
2543  // re-injecting back through producer
2544  cq_s_polymw *= 0.;
2545  }
2546  }
2547  connectionRates[perf][Indices::contiPolymerMWEqIdx] = Base::restrictEval(cq_s_polymw);
2548  }
2549 
2550 
2551 
2552 
2553 
2554 
2555  template<typename TypeTag>
2556  std::optional<double>
2557  StandardWell<TypeTag>::
2558  computeBhpAtThpLimitProd(const WellState& well_state,
2559  const Simulator& ebos_simulator,
2560  const SummaryState& summary_state,
2561  DeferredLogger& deferred_logger) const
2562  {
2563  return computeBhpAtThpLimitProdWithAlq(ebos_simulator,
2564  summary_state,
2565  deferred_logger,
2566  this->getALQ(well_state));
2567  }
2568 
2569  template<typename TypeTag>
2570  std::optional<double>
2571  StandardWell<TypeTag>::
2572  computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
2573  const SummaryState& summary_state,
2574  DeferredLogger& deferred_logger,
2575  double alq_value) const
2576  {
2577  // Make the frates() function.
2578  auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
2579  // Not solving the well equations here, which means we are
2580  // calculating at the current Fg/Fw values of the
2581  // well. This does not matter unless the well is
2582  // crossflowing, and then it is likely still a good
2583  // approximation.
2584  std::vector<double> rates(3);
2585  computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2586  this->adaptRatesForVFP(rates);
2587  return rates;
2588  };
2589 
2590  double max_pressure = 0.0;
2591  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
2592  const int cell_idx = this->well_cells_[perf];
2593  const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
2594  const auto& fs = int_quants.fluidState();
2595  double pressure_cell = this->getPerfCellPressure(fs).value();
2596  max_pressure = std::max(max_pressure, pressure_cell);
2597  }
2598  auto bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(frates,
2599  summary_state,
2600  deferred_logger,
2601  max_pressure,
2602  alq_value);
2603  auto v = frates(*bhpAtLimit);
2604  if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
2605  return bhpAtLimit;
2606 
2607  auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) {
2608  // Solver the well iterations to see if we are
2609  // able to get a solution with an update
2610  // solution
2611  std::vector<double> rates(3);
2612  computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger);
2613  this->adaptRatesForVFP(rates);
2614  return rates;
2615  };
2616 
2617  bhpAtLimit = this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitProdWithAlq(fratesIter,
2618  summary_state,
2619  deferred_logger,
2620  max_pressure,
2621  alq_value);
2622  v = frates(*bhpAtLimit);
2623  if(bhpAtLimit && std::all_of(v.cbegin(), v.cend(), [](double i){ return i <= 0; }))
2624  return bhpAtLimit;
2625 
2626  // we still don't get a valied solution.
2627  return std::nullopt;
2628  }
2629 
2630 
2631 
2632  template<typename TypeTag>
2633  std::optional<double>
2634  StandardWell<TypeTag>::
2635  computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
2636  const SummaryState& summary_state,
2637  DeferredLogger& deferred_logger) const
2638  {
2639  // Make the frates() function.
2640  auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
2641  // Not solving the well equations here, which means we are
2642  // calculating at the current Fg/Fw values of the
2643  // well. This does not matter unless the well is
2644  // crossflowing, and then it is likely still a good
2645  // approximation.
2646  std::vector<double> rates(3);
2647  computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
2648  return rates;
2649  };
2650 
2651  return this->StandardWellGeneric<Scalar>::computeBhpAtThpLimitInj(frates,
2652  summary_state,
2653  deferred_logger);
2654  }
2655 
2656 
2657 
2658 
2659 
2660  template<typename TypeTag>
2661  bool
2662  StandardWell<TypeTag>::
2663  iterateWellEqWithControl(const Simulator& ebosSimulator,
2664  const double dt,
2665  const Well::InjectionControls& inj_controls,
2666  const Well::ProductionControls& prod_controls,
2667  WellState& well_state,
2668  const GroupState& group_state,
2669  DeferredLogger& deferred_logger)
2670  {
2671  const int max_iter = this->param_.max_inner_iter_wells_;
2672  int it = 0;
2673  bool converged;
2674  bool relax_convergence = false;
2675  this->regularize_ = false;
2676  do {
2677  assembleWellEqWithoutIteration(ebosSimulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger);
2678 
2679  if (it > this->param_.strict_inner_iter_wells_) {
2680  relax_convergence = true;
2681  this->regularize_ = true;
2682  }
2683 
2684  auto report = getWellConvergence(well_state, Base::B_avg_, deferred_logger, relax_convergence);
2685 
2686  converged = report.converged();
2687  if (converged) {
2688  break;
2689  }
2690 
2691  ++it;
2692  solveEqAndUpdateWellState(well_state, deferred_logger);
2693 
2694  // TODO: when this function is used for well testing purposes, will need to check the controls, so that we will obtain convergence
2695  // under the most restrictive control. Based on this converged results, we can check whether to re-open the well. Either we refactor
2696  // this function or we use different functions for the well testing purposes.
2697  // We don't allow for switching well controls while computing well potentials and testing wells
2698  // updateWellControl(ebosSimulator, well_state, deferred_logger);
2699  initPrimaryVariablesEvaluation();
2700  } while (it < max_iter);
2701 
2702  return converged;
2703  }
2704 
2705 
2706  template<typename TypeTag>
2707  std::vector<double>
2709  computeCurrentWellRates(const Simulator& ebosSimulator,
2710  DeferredLogger& deferred_logger) const
2711  {
2712  // Calculate the rates that follow from the current primary variables.
2713  std::vector<double> well_q_s(this->num_components_, 0.);
2714  const EvalWell& bhp = this->getBhp();
2715  const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
2716  for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
2717  const int cell_idx = this->well_cells_[perf];
2718  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
2719  std::vector<Scalar> mob(this->num_components_, 0.);
2720  getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
2721  std::vector<Scalar> cq_s(this->num_components_, 0.);
2722  double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
2723  const double Tw = this->well_index_[perf] * trans_mult;
2724  computePerfRateScalar(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
2725  cq_s, deferred_logger);
2726  for (int comp = 0; comp < this->num_components_; ++comp) {
2727  well_q_s[comp] += cq_s[comp];
2728  }
2729  }
2730  const auto& comm = this->parallel_well_info_.communication();
2731  if (comm.size() > 1)
2732  {
2733  comm.sum(well_q_s.data(), well_q_s.size());
2734  }
2735  return well_q_s;
2736  }
2737 
2738 
2739 
2740 
2741 
2742  template <typename TypeTag>
2743  void
2745  computeConnLevelProdInd(const typename StandardWell<TypeTag>::FluidState& fs,
2746  const std::function<double(const double)>& connPICalc,
2747  const std::vector<EvalWell>& mobility,
2748  double* connPI) const
2749  {
2750  const auto& pu = this->phaseUsage();
2751  const int np = this->number_of_phases_;
2752  for (int p = 0; p < np; ++p) {
2753  // Note: E100's notion of PI value phase mobility includes
2754  // the reciprocal FVF.
2755  const auto connMob =
2756  mobility[ this->flowPhaseToEbosCompIdx(p) ].value()
2757  * fs.invB(this->flowPhaseToEbosPhaseIdx(p)).value();
2758 
2759  connPI[p] = connPICalc(connMob);
2760  }
2761 
2762  if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
2763  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
2764  {
2765  const auto io = pu.phase_pos[Oil];
2766  const auto ig = pu.phase_pos[Gas];
2767 
2768  const auto vapoil = connPI[ig] * fs.Rv().value();
2769  const auto disgas = connPI[io] * fs.Rs().value();
2770 
2771  connPI[io] += vapoil;
2772  connPI[ig] += disgas;
2773  }
2774  }
2775 
2776 
2777 
2778 
2779 
2780  template <typename TypeTag>
2781  void
2782  StandardWell<TypeTag>::
2783  computeConnLevelInjInd(const typename StandardWell<TypeTag>::FluidState& fs,
2784  const Phase preferred_phase,
2785  const std::function<double(const double)>& connIICalc,
2786  const std::vector<EvalWell>& mobility,
2787  double* connII,
2788  DeferredLogger& deferred_logger) const
2789  {
2790  // Assumes single phase injection
2791  const auto& pu = this->phaseUsage();
2792 
2793  auto phase_pos = 0;
2794  if (preferred_phase == Phase::GAS) {
2795  phase_pos = pu.phase_pos[Gas];
2796  }
2797  else if (preferred_phase == Phase::OIL) {
2798  phase_pos = pu.phase_pos[Oil];
2799  }
2800  else if (preferred_phase == Phase::WATER) {
2801  phase_pos = pu.phase_pos[Water];
2802  }
2803  else {
2804  OPM_DEFLOG_THROW(NotImplemented,
2805  "Unsupported Injector Type ("
2806  << static_cast<int>(preferred_phase)
2807  << ") for well " << this->name()
2808  << " during connection I.I. calculation",
2809  deferred_logger);
2810  }
2811 
2812  const auto zero = EvalWell { this->numWellEq_ + Indices::numEq, 0.0 };
2813  const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
2814  connII[phase_pos] = connIICalc(mt.value() * fs.invB(this->flowPhaseToEbosPhaseIdx(phase_pos)).value());
2815  }
2816 } // 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: StandardWell.hpp:63
Collect per-connection static information to enable calculating connection-level or well-level produc...
Definition: WellProdIndexCalculator.hpp:36
double connectionProdIndStandard(const std::size_t connIdx, const double connMobility) const
Compute connection-level steady-state productivity index value using dynamic phase mobility.
Definition: WellProdIndexCalculator.cpp:110
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