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