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