All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
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 
23 namespace Opm
24 {
25  template<typename TypeTag>
26  StandardWell<TypeTag>::
27  StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
28  : Base(well, time_step, wells, param)
29  , perf_densities_(number_of_perforations_)
30  , perf_pressure_diffs_(number_of_perforations_)
31  , primary_variables_(numWellEq, 0.0)
32  , primary_variables_evaluation_(numWellEq) // the number of the primary variables
33  , F0_(numWellEq)
34  {
35  duneB_.setBuildMode( OffDiagMatWell::row_wise );
36  duneC_.setBuildMode( OffDiagMatWell::row_wise );
37  invDuneD_.setBuildMode( DiagMatWell::row_wise );
38  }
39 
40 
41 
42 
43 
44  template<typename TypeTag>
45  void
46  StandardWell<TypeTag>::
47  init(const PhaseUsage* phase_usage_arg,
48  const std::vector<bool>* active_arg,
49  const std::vector<double>& depth_arg,
50  const double gravity_arg,
51  const int num_cells)
52  {
53  Base::init(phase_usage_arg, active_arg,
54  depth_arg, gravity_arg, num_cells);
55 
56  perf_depth_.resize(number_of_perforations_, 0.);
57  for (int perf = 0; perf < number_of_perforations_; ++perf) {
58  const int cell_idx = well_cells_[perf];
59  perf_depth_[perf] = depth_arg[cell_idx];
60  }
61 
62  // setup sparsity pattern for the matrices
63  //[A C^T [x = [ res
64  // B D] x_well] res_well]
65  // set the size of the matrices
66  invDuneD_.setSize(1, 1, 1);
67  duneB_.setSize(1, num_cells, number_of_perforations_);
68  duneC_.setSize(1, num_cells, number_of_perforations_);
69 
70  for (auto row=invDuneD_.createbegin(), end = invDuneD_.createend(); row!=end; ++row) {
71  // Add nonzeros for diagonal
72  row.insert(row.index());
73  }
74 
75  for (auto row = duneB_.createbegin(), end = duneB_.createend(); row!=end; ++row) {
76  for (int perf = 0 ; perf < number_of_perforations_; ++perf) {
77  const int cell_idx = well_cells_[perf];
78  row.insert(cell_idx);
79  }
80  }
81 
82  // make the C^T matrix
83  for (auto row = duneC_.createbegin(), end = duneC_.createend(); row!=end; ++row) {
84  for (int perf = 0; perf < number_of_perforations_; ++perf) {
85  const int cell_idx = well_cells_[perf];
86  row.insert(cell_idx);
87  }
88  }
89 
90  resWell_.resize(1);
91 
92  // resize temporary class variables
93  Bx_.resize( duneB_.N() );
94  invDrw_.resize( invDuneD_.N() );
95  }
96 
97 
98 
99 
100 
101  template<typename TypeTag>
102  void StandardWell<TypeTag>::
103  initPrimaryVariablesEvaluation() const
104  {
105  // TODO: using numComp here is only to make the 2p + dummy phase work
106  // TODO: in theory, we should use numWellEq here.
107  // for (int eqIdx = 0; eqIdx < numWellEq; ++eqIdx) {
108  for (int eqIdx = 0; eqIdx < numComponents(); ++eqIdx) {
109  assert( (size_t)eqIdx < primary_variables_.size() );
110 
111  primary_variables_evaluation_[eqIdx] = 0.0;
112  primary_variables_evaluation_[eqIdx].setValue(primary_variables_[eqIdx]);
113  primary_variables_evaluation_[eqIdx].setDerivative(numEq + eqIdx, 1.0);
114  }
115  }
116 
117 
118 
119 
120 
121  template<typename TypeTag>
122  typename StandardWell<TypeTag>::EvalWell
123  StandardWell<TypeTag>::
124  getBhp() const
125  {
126  const WellControls* wc = well_controls_;
127  if (well_controls_get_current_type(wc) == BHP) {
128  EvalWell bhp = 0.0;
129  const double target_rate = well_controls_get_current_target(wc);
130  bhp.setValue(target_rate);
131  return bhp;
132  } else if (well_controls_get_current_type(wc) == THP) {
133  const int control = well_controls_get_current(wc);
134 
135  const Opm::PhaseUsage& pu = phaseUsage();
136  std::vector<EvalWell> rates(3, 0.0);
137  if (active()[ Water ]) {
138  rates[ Water ]= getQs(pu.phase_pos[ Water]);
139  }
140  if (active()[ Oil ]) {
141  rates[ Oil ] = getQs(pu.phase_pos[ Oil ]);
142  }
143  if (active()[ Gas ]) {
144  rates[ Gas ] = getQs(pu.phase_pos[ Gas ]);
145  }
146  return calculateBhpFromThp(rates, control);
147  }
148 
149  return primary_variables_evaluation_[XvarWell];
150  }
151 
152 
153 
154 
155 
156  template<typename TypeTag>
157  typename StandardWell<TypeTag>::EvalWell
158  StandardWell<TypeTag>::
159  getQs(const int comp_idx) const // TODO: phase or component?
160  {
161  EvalWell qs = 0.0;
162 
163  const WellControls* wc = well_controls_;
164  const int np = number_of_phases_;
165  const double target_rate = well_controls_get_current_target(wc);
166 
167  assert(comp_idx < numComponents());
168  const auto pu = phaseUsage();
169 
170  // TODO: the formulation for the injectors decides it only work with single phase
171  // surface rate injection control. Improvement will be required.
172  if (well_type_ == INJECTOR) {
173  if (has_solvent) {
174  // TODO: investigate whether the use of the comp_frac is justified.
175  // The usage of the comp_frac is not correct, which should be changed later.
176  double comp_frac = 0.0;
177  if (has_solvent && comp_idx == contiSolventEqIdx) { // solvent
178  comp_frac = comp_frac_[pu.phase_pos[ Gas ]] * wsolvent();
179  } else if (comp_idx == pu.phase_pos[ Gas ]) {
180  comp_frac = comp_frac_[comp_idx] * (1.0 - wsolvent());
181  } else {
182  comp_frac = comp_frac_[comp_idx];
183  }
184  if (comp_frac == 0.0) {
185  return qs; //zero
186  }
187 
188  if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
189  return comp_frac * primary_variables_evaluation_[XvarWell];
190  }
191 
192  qs.setValue(comp_frac * target_rate);
193  return qs;
194  }
195 
196  const double comp_frac = comp_frac_[comp_idx];
197  if (comp_frac == 0.0) {
198  return qs;
199  }
200 
201  if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
202  return primary_variables_evaluation_[XvarWell];
203  }
204  qs.setValue(target_rate);
205  return qs;
206  }
207 
208  // Producers
209  if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
210  return primary_variables_evaluation_[XvarWell] * wellVolumeFractionScaled(comp_idx);
211  }
212 
213  if (well_controls_get_current_type(wc) == SURFACE_RATE) {
214  // checking how many phases are included in the rate control
215  // to decide wheter it is a single phase rate control or not
216  const double* distr = well_controls_get_current_distr(wc);
217  int num_phases_under_rate_control = 0;
218  for (int phase = 0; phase < np; ++phase) {
219  if (distr[phase] > 0.0) {
220  num_phases_under_rate_control += 1;
221  }
222  }
223 
224  // there should be at least one phase involved
225  assert(num_phases_under_rate_control > 0);
226 
227  // when it is a single phase rate limit
228  if (num_phases_under_rate_control == 1) {
229 
230  // looking for the phase under control
231  int phase_under_control = -1;
232  for (int phase = 0; phase < np; ++phase) {
233  if (distr[phase] > 0.0) {
234  phase_under_control = phase;
235  break;
236  }
237  }
238 
239  assert(phase_under_control >= 0);
240 
241  EvalWell wellVolumeFractionScaledPhaseUnderControl = wellVolumeFractionScaled(phase_under_control);
242  if (has_solvent && phase_under_control == Gas) {
243  // for GRAT controlled wells solvent is included in the target
244  wellVolumeFractionScaledPhaseUnderControl += wellVolumeFractionScaled(contiSolventEqIdx);
245  }
246 
247  if (comp_idx == phase_under_control) {
248  if (has_solvent && phase_under_control == Gas) {
249  qs.setValue(target_rate * wellVolumeFractionScaled(Gas).value() / wellVolumeFractionScaledPhaseUnderControl.value() );
250  return qs;
251  }
252  qs.setValue(target_rate);
253  return qs;
254  }
255 
256  // TODO: not sure why the single phase under control will have near zero fraction
257  const double eps = 1e-6;
258  if (wellVolumeFractionScaledPhaseUnderControl < eps) {
259  return qs;
260  }
261  return (target_rate * wellVolumeFractionScaled(comp_idx) / wellVolumeFractionScaledPhaseUnderControl);
262  }
263 
264  // when it is a combined two phase rate limit, such like LRAT
265  // we neec to calculate the rate for the certain phase
266  if (num_phases_under_rate_control == 2) {
267  EvalWell combined_volume_fraction = 0.;
268  for (int p = 0; p < np; ++p) {
269  if (distr[p] == 1.0) {
270  combined_volume_fraction += wellVolumeFractionScaled(p);
271  }
272  }
273  return (target_rate * wellVolumeFractionScaled(comp_idx) / combined_volume_fraction);
274  }
275 
276  // TODO: three phase surface rate control is not tested yet
277  if (num_phases_under_rate_control == 3) {
278  return target_rate * wellSurfaceVolumeFraction(comp_idx);
279  }
280  } else if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
281  // ReservoirRate
282  return target_rate * wellVolumeFractionScaled(comp_idx);
283  } else {
284  OPM_THROW(std::logic_error, "Unknown control type for well " << name());
285  }
286 
287  // avoid warning of condition reaches end of non-void function
288  return qs;
289  }
290 
291 
292 
293 
294 
295 
296  template<typename TypeTag>
297  typename StandardWell<TypeTag>::EvalWell
298  StandardWell<TypeTag>::
299  wellVolumeFractionScaled(const int compIdx) const
300  {
301 
302  const double scal = scalingFactor(compIdx);
303  if (scal > 0)
304  return wellVolumeFraction(compIdx) / scal;
305 
306  // the scaling factor may be zero for RESV controlled wells.
307  return wellVolumeFraction(compIdx);
308  }
309 
310 
311 
312 
313 
314  template<typename TypeTag>
315  typename StandardWell<TypeTag>::EvalWell
316  StandardWell<TypeTag>::
317  wellVolumeFraction(const int compIdx) const
318  {
319  const auto& pu = phaseUsage();
320  if (active()[Water] && compIdx == pu.phase_pos[Water]) {
321  return primary_variables_evaluation_[WFrac];
322  }
323 
324  if (active()[Gas] && compIdx == pu.phase_pos[Gas]) {
325  return primary_variables_evaluation_[GFrac];
326  }
327 
328  if (has_solvent && compIdx == contiSolventEqIdx) {
329  return primary_variables_evaluation_[SFrac];
330  }
331 
332  // Oil fraction
333  EvalWell well_fraction = 1.0;
334  if (active()[Water]) {
335  well_fraction -= primary_variables_evaluation_[WFrac];
336  }
337 
338  if (active()[Gas]) {
339  well_fraction -= primary_variables_evaluation_[GFrac];
340  }
341  if (has_solvent) {
342  well_fraction -= primary_variables_evaluation_[SFrac];
343  }
344  return well_fraction;
345  }
346 
347 
348 
349 
350 
351  template<typename TypeTag>
352  typename StandardWell<TypeTag>::EvalWell
353  StandardWell<TypeTag>::
354  wellSurfaceVolumeFraction(const int compIdx) const
355  {
356  EvalWell sum_volume_fraction_scaled = 0.;
357  const int numComp = numComponents();
358  for (int idx = 0; idx < numComp; ++idx) {
359  sum_volume_fraction_scaled += wellVolumeFractionScaled(idx);
360  }
361 
362  assert(sum_volume_fraction_scaled.value() != 0.);
363 
364  return wellVolumeFractionScaled(compIdx) / sum_volume_fraction_scaled;
365  }
366 
367 
368 
369 
370 
371  template<typename TypeTag>
372  typename StandardWell<TypeTag>::EvalWell
373  StandardWell<TypeTag>::
374  extendEval(const Eval& in) const
375  {
376  EvalWell out = 0.0;
377  out.setValue(in.value());
378  for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
379  out.setDerivative(eqIdx, in.derivative(eqIdx));
380  }
381  return out;
382  }
383 
384 
385 
386 
387 
388  template<typename TypeTag>
389  void
390  StandardWell<TypeTag>::
391  computePerfRate(const IntensiveQuantities& intQuants,
392  const std::vector<EvalWell>& mob_perfcells_dense,
393  const double Tw, const EvalWell& bhp, const double& cdp,
394  const bool& allow_cf, std::vector<EvalWell>& cq_s) const
395  {
396  const Opm::PhaseUsage& pu = phaseUsage();
397  const int np = number_of_phases_;
398  const int numComp = numComponents();
399  std::vector<EvalWell> cmix_s(numComp,0.0);
400  for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
401  cmix_s[componentIdx] = wellSurfaceVolumeFraction(componentIdx);
402  }
403  const auto& fs = intQuants.fluidState();
404 
405  const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
406  const EvalWell rs = extendEval(fs.Rs());
407  const EvalWell rv = extendEval(fs.Rv());
408  std::vector<EvalWell> b_perfcells_dense(numComp, 0.0);
409  for (int phase = 0; phase < np; ++phase) {
410  const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
411  b_perfcells_dense[phase] = extendEval(fs.invB(ebosPhaseIdx));
412  }
413  if (has_solvent) {
414  b_perfcells_dense[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
415  }
416 
417  // Pressure drawdown (also used to determine direction of flow)
418  const EvalWell well_pressure = bhp + cdp;
419  const EvalWell drawdown = pressure - well_pressure;
420 
421  // producing perforations
422  if ( drawdown.value() > 0 ) {
423  //Do nothing if crossflow is not allowed
424  if (!allow_cf && well_type_ == INJECTOR) {
425  return;
426  }
427 
428  // compute component volumetric rates at standard conditions
429  for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
430  const EvalWell cq_p = - Tw * (mob_perfcells_dense[componentIdx] * drawdown);
431  cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
432  }
433 
434  if (active()[Oil] && active()[Gas]) {
435  const int oilpos = pu.phase_pos[Oil];
436  const int gaspos = pu.phase_pos[Gas];
437  const EvalWell cq_sOil = cq_s[oilpos];
438  const EvalWell cq_sGas = cq_s[gaspos];
439  cq_s[gaspos] += rs * cq_sOil;
440  cq_s[oilpos] += rv * cq_sGas;
441  }
442 
443  } else {
444  //Do nothing if crossflow is not allowed
445  if (!allow_cf && well_type_ == PRODUCER) {
446  return;
447  }
448 
449  // Using total mobilities
450  EvalWell total_mob_dense = mob_perfcells_dense[0];
451  for (int componentIdx = 1; componentIdx < numComp; ++componentIdx) {
452  total_mob_dense += mob_perfcells_dense[componentIdx];
453  }
454 
455  // injection perforations total volume rates
456  const EvalWell cqt_i = - Tw * (total_mob_dense * drawdown);
457 
458  // compute volume ratio between connection at standard conditions
459  EvalWell volumeRatio = 0.0;
460  if (active()[Water]) {
461  const int watpos = pu.phase_pos[Water];
462  volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos];
463  }
464 
465  if (has_solvent) {
466  volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx];
467  }
468 
469  if (active()[Oil] && active()[Gas]) {
470  const int oilpos = pu.phase_pos[Oil];
471  const int gaspos = pu.phase_pos[Gas];
472 
473  // Incorporate RS/RV factors if both oil and gas active
474  const EvalWell d = 1.0 - rv * rs;
475 
476  if (d.value() == 0.0) {
477  OPM_THROW(Opm::NumericalProblem, "Zero d value obtained for well " << name() << " during flux calcuation"
478  << " with rs " << rs << " and rv " << rv);
479  }
480 
481  const EvalWell tmp_oil = (cmix_s[oilpos] - rv * cmix_s[gaspos]) / d;
482  //std::cout << "tmp_oil " <<tmp_oil << std::endl;
483  volumeRatio += tmp_oil / b_perfcells_dense[oilpos];
484 
485  const EvalWell tmp_gas = (cmix_s[gaspos] - rs * cmix_s[oilpos]) / d;
486  //std::cout << "tmp_gas " <<tmp_gas << std::endl;
487  volumeRatio += tmp_gas / b_perfcells_dense[gaspos];
488  }
489  else {
490  if (active()[Oil]) {
491  const int oilpos = pu.phase_pos[Oil];
492  volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos];
493  }
494  if (active()[Gas]) {
495  const int gaspos = pu.phase_pos[Gas];
496  volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos];
497  }
498  }
499 
500  // injecting connections total volumerates at standard conditions
501  EvalWell cqt_is = cqt_i/volumeRatio;
502  //std::cout << "volrat " << volumeRatio << " " << volrat_perf_[perf] << std::endl;
503  for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
504  cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; // * b_perfcells_dense[phase];
505  }
506  }
507  }
508 
509 
510 
511 
512 
513  template<typename TypeTag>
514  void
515  StandardWell<TypeTag>::
516  assembleWellEq(Simulator& ebosSimulator,
517  const double dt,
518  WellState& well_state,
519  bool only_wells)
520  {
521  const int numComp = numComponents();
522  const int np = number_of_phases_;
523 
524  // clear all entries
525  if (!only_wells) {
526  duneB_ = 0.0;
527  duneC_ = 0.0;
528  }
529  invDuneD_ = 0.0;
530  resWell_ = 0.0;
531 
532  auto& ebosJac = ebosSimulator.model().linearizer().matrix();
533  auto& ebosResid = ebosSimulator.model().linearizer().residual();
534 
535  // TODO: it probably can be static member for StandardWell
536  const double volume = 0.002831684659200; // 0.1 cu ft;
537 
538  const bool allow_cf = crossFlowAllowed(ebosSimulator);
539 
540  const EvalWell& bhp = getBhp();
541 
542  for (int perf = 0; perf < number_of_perforations_; ++perf) {
543 
544  const int cell_idx = well_cells_[perf];
545  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
546  std::vector<EvalWell> cq_s(numComp,0.0);
547  std::vector<EvalWell> mob(numComp, 0.0);
548  getMobility(ebosSimulator, perf, mob);
549  computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
550 
551  for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
552  // the cq_s entering mass balance equations need to consider the efficiency factors.
553  const EvalWell cq_s_effective = cq_s[componentIdx] * well_efficiency_factor_;
554 
555  if (!only_wells) {
556  // subtract sum of component fluxes in the reservoir equation.
557  // need to consider the efficiency factor
558  ebosResid[cell_idx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.value();
559  }
560 
561  // subtract sum of phase fluxes in the well equations.
562  resWell_[0][componentIdx] -= cq_s_effective.value();
563 
564  // assemble the jacobians
565  for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
566  if (!only_wells) {
567  // also need to consider the efficiency factor when manipulating the jacobians.
568  duneC_[0][cell_idx][pvIdx][flowPhaseToEbosCompIdx(componentIdx)] -= cq_s_effective.derivative(pvIdx+numEq); // intput in transformed matrix
569  }
570  invDuneD_[0][0][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx+numEq);
571  }
572 
573  for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
574  if (!only_wells) {
575  // also need to consider the efficiency factor when manipulating the jacobians.
576  ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(componentIdx)][pvIdx] -= cq_s_effective.derivative(pvIdx);
577  duneB_[0][cell_idx][componentIdx][pvIdx] -= cq_s_effective.derivative(pvIdx);
578  }
579  }
580 
581  // Store the perforation phase flux for later usage.
582  if (has_solvent && componentIdx == contiSolventEqIdx) {// if (flowPhaseToEbosCompIdx(componentIdx) == Solvent)
583  well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value();
584  } else {
585  well_state.perfPhaseRates()[(first_perf_ + perf) * np + componentIdx] = cq_s[componentIdx].value();
586  }
587  }
588 
589  if (has_polymer) {
590  EvalWell cq_s_poly = cq_s[Water];
591  if (well_type_ == INJECTOR) {
592  cq_s_poly *= wpolymer();
593  } else {
594  cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
595  }
596  if (!only_wells) {
597  // TODO: we need to consider the efficiency here.
598  for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
599  ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
600  }
601  ebosResid[cell_idx][contiPolymerEqIdx] -= cq_s_poly.value();
602  }
603  }
604 
605  // Store the perforation pressure for later usage.
606  well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
607  }
608 
609  // add vol * dF/dt + Q to the well equations;
610  for (int componentIdx = 0; componentIdx < numComp; ++componentIdx) {
611  EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;
612  resWell_loc += getQs(componentIdx) * well_efficiency_factor_;
613  for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
614  invDuneD_[0][0][componentIdx][pvIdx] += resWell_loc.derivative(pvIdx+numEq);
615  }
616  resWell_[0][componentIdx] += resWell_loc.value();
617  }
618 
619  // do the local inversion of D.
620  invDuneD_[0][0].invert();
621  }
622 
623 
624 
625 
626 
627  template<typename TypeTag>
628  bool
629  StandardWell<TypeTag>::
630  crossFlowAllowed(const Simulator& ebosSimulator) const
631  {
632  if (getAllowCrossFlow()) {
633  return true;
634  }
635 
636  // TODO: investigate the justification of the following situation
637 
638  // check for special case where all perforations have cross flow
639  // then the wells must allow for cross flow
640  for (int perf = 0; perf < number_of_perforations_; ++perf) {
641  const int cell_idx = well_cells_[perf];
642  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
643  const auto& fs = intQuants.fluidState();
644  EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
645  EvalWell bhp = getBhp();
646 
647  // Pressure drawdown (also used to determine direction of flow)
648  EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
649  EvalWell drawdown = pressure - well_pressure;
650 
651  if (drawdown.value() < 0 && well_type_ == INJECTOR) {
652  return false;
653  }
654 
655  if (drawdown.value() > 0 && well_type_ == PRODUCER) {
656  return false;
657  }
658  }
659  return true;
660  }
661 
662 
663 
664 
665 
666  template<typename TypeTag>
667  void
668  StandardWell<TypeTag>::
669  getMobility(const Simulator& ebosSimulator,
670  const int perf,
671  std::vector<EvalWell>& mob) const
672  {
673  const int np = number_of_phases_;
674  const int cell_idx = well_cells_[perf];
675  assert (int(mob.size()) == numComponents());
676  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
677  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
678 
679  // either use mobility of the perforation cell or calcualte its own
680  // based on passing the saturation table index
681  const int satid = saturation_table_number_[perf] - 1;
682  const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx);
683  if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell
684 
685  for (int phase = 0; phase < np; ++phase) {
686  int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
687  mob[phase] = extendEval(intQuants.mobility(ebosPhaseIdx));
688  }
689  if (has_solvent) {
690  mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility());
691  }
692  } else {
693 
694  const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
695  Eval relativePerms[3] = { 0.0, 0.0, 0.0 };
696  MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState());
697 
698  // reset the satnumvalue back to original
699  materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx);
700 
701  // compute the mobility
702  for (int phase = 0; phase < np; ++phase) {
703  int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phase);
704  mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
705  }
706 
707  // this may not work if viscosity and relperms has been modified?
708  if (has_solvent) {
709  OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent");
710  }
711  }
712 
713  // modify the water mobility if polymer is present
714  if (has_polymer) {
715  // assume fully mixture for wells.
716  EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
717 
718  if (well_type_ == INJECTOR) {
719  const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex());
720  mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) );
721  }
722 
723  if (PolymerModule::hasPlyshlog()) {
724  // compute the well water velocity with out shear effects.
725  const int numComp = numComponents();
726  const bool allow_cf = crossFlowAllowed(ebosSimulator);
727  const EvalWell& bhp = getBhp();
728  std::vector<EvalWell> cq_s(numComp,0.0);
729  computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
730  // TODO: make area a member
731  double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf];
732  const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
733  const auto& scaledDrainageInfo =
734  materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx);
735  const Scalar& Swcr = scaledDrainageInfo.Swcr;
736  const EvalWell poro = extendEval(intQuants.porosity());
737  const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water)));
738  // guard against zero porosity and no water
739  const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12);
740  EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water)));
741 
742  if (PolymerModule::hasShrate()) {
743  // TODO Use the same conversion as for the reservoar equations.
744  // Need the "permeability" of the well?
745  // For now use the same formula as in legacy.
746  waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / bore_diameters_[perf];
747  }
748  EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
749  EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration,
750  intQuants.pvtRegionIndex(),
751  waterVelocity);
752 
753  // modify the mobility with the shear factor and recompute the well fluxes.
754  mob[ Water ] /= shearFactor;
755  }
756  }
757  }
758 
759 
760 
761 
762 
763  template<typename TypeTag>
764  void
765  StandardWell<TypeTag>::
766  updateWellState(const BVectorWell& dwells,
767  WellState& well_state) const
768  {
769  const int np = number_of_phases_;
770  const double dBHPLimit = param_.dbhp_max_rel_;
771  const double dFLimit = param_.dwell_fraction_max_;
772  const auto pu = phaseUsage();
773 
774  const std::vector<double> xvar_well_old = primary_variables_;
775 
776  // update the second and third well variable (The flux fractions)
777  std::vector<double> F(np,0.0);
778  if (active()[ Water ]) {
779  const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
780  const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac]),dFLimit);
781  primary_variables_[WFrac] = xvar_well_old[WFrac] - dx2_limited;
782  }
783 
784  if (active()[ Gas ]) {
785  const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
786  const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac]),dFLimit);
787  primary_variables_[GFrac] = xvar_well_old[GFrac] - dx3_limited;
788  }
789 
790  if (has_solvent) {
791  const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
792  const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]),dFLimit);
793  primary_variables_[SFrac] = xvar_well_old[SFrac] - dx4_limited;
794  }
795 
796  assert(active()[ Oil ]);
797  F[pu.phase_pos[Oil]] = 1.0;
798 
799  if (active()[ Water ]) {
800  F[pu.phase_pos[Water]] = primary_variables_[WFrac];
801  F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]];
802  }
803 
804  if (active()[ Gas ]) {
805  F[pu.phase_pos[Gas]] = primary_variables_[GFrac];
806  F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Gas]];
807  }
808 
809  double F_solvent = 0.0;
810  if (has_solvent) {
811  F_solvent = primary_variables_[SFrac];
812  F[pu.phase_pos[Oil]] -= F_solvent;
813  }
814 
815  if (active()[ Water ]) {
816  if (F[Water] < 0.0) {
817  if (active()[ Gas ]) {
818  F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Water]]);
819  }
820  if (has_solvent) {
821  F_solvent /= (1.0 - F[pu.phase_pos[Water]]);
822  }
823  F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Water]]);
824  F[pu.phase_pos[Water]] = 0.0;
825  }
826  }
827 
828  if (active()[ Gas ]) {
829  if (F[pu.phase_pos[Gas]] < 0.0) {
830  if (active()[ Water ]) {
831  F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Gas]]);
832  }
833  if (has_solvent) {
834  F_solvent /= (1.0 - F[pu.phase_pos[Gas]]);
835  }
836  F[pu.phase_pos[Oil]] /= (1.0 - F[pu.phase_pos[Gas]]);
837  F[pu.phase_pos[Gas]] = 0.0;
838  }
839  }
840 
841  if (F[pu.phase_pos[Oil]] < 0.0) {
842  if (active()[ Water ]) {
843  F[pu.phase_pos[Water]] /= (1.0 - F[pu.phase_pos[Oil]]);
844  }
845  if (active()[ Gas ]) {
846  F[pu.phase_pos[Gas]] /= (1.0 - F[pu.phase_pos[Oil]]);
847  }
848  if (has_solvent) {
849  F_solvent /= (1.0 - F[pu.phase_pos[Oil]]);
850  }
851  F[pu.phase_pos[Oil]] = 0.0;
852  }
853 
854  if (active()[ Water ]) {
855  primary_variables_[WFrac] = F[pu.phase_pos[Water]];
856  }
857  if (active()[ Gas ]) {
858  primary_variables_[GFrac] = F[pu.phase_pos[Gas]];
859  }
860  if(has_solvent) {
861  primary_variables_[SFrac] = F_solvent;
862  }
863 
864  // F_solvent is added to F_gas. This means that well_rate[Gas] also contains solvent.
865  // More testing is needed to make sure this is correct for well groups and THP.
866  if (has_solvent){
867  F[pu.phase_pos[Gas]] += F_solvent;
868  }
869 
870  // The interpretation of the first well variable depends on the well control
871  const WellControls* wc = well_controls_;
872 
873  // TODO: we should only maintain one current control either from the well_state or from well_controls struct.
874  // Either one can be more favored depending on the final strategy for the initilzation of the well control
875  const int current = well_state.currentControls()[index_of_well_];
876  const double target_rate = well_controls_iget_target(wc, current);
877 
878  for (int p = 0; p < np; ++p) {
879  const double scal = scalingFactor(p);
880  if (scal > 0) {
881  F[p] /= scal ;
882  } else {
883  F[p] = 0.;
884  }
885  }
886 
887  switch (well_controls_iget_type(wc, current)) {
888  case THP: // The BHP and THP both uses the total rate as first well variable.
889  case BHP:
890  {
891  primary_variables_[XvarWell] = xvar_well_old[XvarWell] - dwells[0][XvarWell];
892 
893  switch (well_type_) {
894  case INJECTOR:
895  for (int p = 0; p < np; ++p) {
896  const double comp_frac = comp_frac_[p];
897  well_state.wellRates()[index_of_well_ * np + p] = comp_frac * primary_variables_[XvarWell];
898  }
899  break;
900  case PRODUCER:
901  for (int p = 0; p < np; ++p) {
902  well_state.wellRates()[index_of_well_ * np + p] = primary_variables_[XvarWell] * F[p];
903  }
904  break;
905  }
906 
907  if (well_controls_iget_type(wc, current) == THP) {
908 
909  // Calculate bhp from thp control and well rates
910  std::vector<double> rates(3, 0.0); // the vfp related only supports three phases for the moment
911 
912  const Opm::PhaseUsage& pu = phaseUsage();
913  if (active()[ Water ]) {
914  rates[ Water ] = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ];
915  }
916  if (active()[ Oil ]) {
917  rates[ Oil ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
918  }
919  if (active()[ Gas ]) {
920  rates[ Gas ]= well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Gas ] ];
921  }
922 
923  well_state.bhp()[index_of_well_] = calculateBhpFromThp(rates, current);
924  }
925  }
926  break;
927  case SURFACE_RATE: // Both rate controls use bhp as first well variable
928  case RESERVOIR_RATE:
929  {
930  const int sign1 = dwells[0][XvarWell] > 0 ? 1: -1;
931  const double dx1_limited = sign1 * std::min(std::abs(dwells[0][XvarWell]),std::abs(xvar_well_old[XvarWell])*dBHPLimit);
932  primary_variables_[XvarWell] = std::max(xvar_well_old[XvarWell] - dx1_limited,1e5);
933  well_state.bhp()[index_of_well_] = primary_variables_[XvarWell];
934 
935  if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
936  if (well_type_ == PRODUCER) {
937 
938  const double* distr = well_controls_iget_distr(wc, current);
939 
940  double F_target = 0.0;
941  for (int p = 0; p < np; ++p) {
942  F_target += distr[p] * F[p];
943  }
944  for (int p = 0; p < np; ++p) {
945  well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate / F_target;
946  }
947  } else {
948 
949  for (int p = 0; p < np; ++p) {
950  well_state.wellRates()[index_of_well_ * np + p] = comp_frac_[p] * target_rate;
951  }
952  }
953  } else { // RESERVOIR_RATE
954  for (int p = 0; p < np; ++p) {
955  well_state.wellRates()[np * index_of_well_ + p] = F[p] * target_rate;
956  }
957  }
958  }
959  break;
960  } // end of switch (well_controls_iget_type(wc, current))
961 
962  // for the wells having a THP constaint, we should update their thp value
963  // If it is under THP control, it will be set to be the target value. Otherwise,
964  // the thp value will be calculated based on the bhp value, assuming the bhp value is correctly calculated.
965  const int nwc = well_controls_get_num(wc);
966  // Looping over all controls until we find a THP constraint
967  int ctrl_index = 0;
968  for ( ; ctrl_index < nwc; ++ctrl_index) {
969  if (well_controls_iget_type(wc, ctrl_index) == THP) {
970  // the current control
971  const int current = well_state.currentControls()[index_of_well_];
972  // If under THP control at the moment
973  if (current == ctrl_index) {
974  const double thp_target = well_controls_iget_target(wc, current);
975  well_state.thp()[index_of_well_] = thp_target;
976  } else { // otherwise we calculate the thp from the bhp value
977  const Opm::PhaseUsage& pu = phaseUsage();
978  std::vector<double> rates(3, 0.0);
979 
980  if (active()[ Water ]) {
981  rates[ Water ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Water ] ];
982  }
983  if (active()[ Oil ]) {
984  rates[ Oil ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Oil ] ];
985  }
986  if (active()[ Gas ]) {
987  rates[ Gas ] = well_state.wellRates()[index_of_well_*np + pu.phase_pos[ Gas ] ];
988  }
989 
990  const double bhp = well_state.bhp()[index_of_well_];
991 
992  well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp);
993  }
994 
995  // the THP control is found, we leave the loop now
996  break;
997  }
998  } // end of for loop for seaching THP constraints
999 
1000  // no THP constraint found
1001  if (ctrl_index == nwc) { // not finding a THP contstraints
1002  well_state.thp()[index_of_well_] = 0.0;
1003  }
1004  }
1005 
1006 
1007 
1008 
1009 
1010  template<typename TypeTag>
1011  void
1013  updateWellStateWithTarget(const int current,
1014  WellState& xw) const
1015  {
1016  // number of phases
1017  const int np = number_of_phases_;
1018  const int well_index = index_of_well_;
1019  const WellControls* wc = well_controls_;
1020  // Updating well state and primary variables.
1021  // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
1022  const double target = well_controls_iget_target(wc, current);
1023  const double* distr = well_controls_iget_distr(wc, current);
1024  switch (well_controls_iget_type(wc, current)) {
1025  case BHP:
1026  xw.bhp()[well_index] = target;
1027  // TODO: similar to the way below to handle THP
1028  // we should not something related to thp here when there is thp constraint
1029  break;
1030 
1031  case THP: {
1032  xw.thp()[well_index] = target;
1033 
1034  const Opm::PhaseUsage& pu = phaseUsage();
1035  std::vector<double> rates(3, 0.0);
1036  if (active()[ Water ]) {
1037  rates[ Water ] = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
1038  }
1039  if (active()[ Oil ]) {
1040  rates[ Oil ] = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
1041  }
1042  if (active()[ Gas ]) {
1043  rates[ Gas ] = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
1044  }
1045 
1046  xw.bhp()[well_index] = calculateBhpFromThp(rates, current);
1047  break;
1048  }
1049 
1050  case RESERVOIR_RATE: // intentional fall-through
1051  case SURFACE_RATE:
1052  // checking the number of the phases under control
1053  int numPhasesWithTargetsUnderThisControl = 0;
1054  for (int phase = 0; phase < np; ++phase) {
1055  if (distr[phase] > 0.0) {
1056  numPhasesWithTargetsUnderThisControl += 1;
1057  }
1058  }
1059 
1060  assert(numPhasesWithTargetsUnderThisControl > 0);
1061 
1062  if (well_type_ == INJECTOR) {
1063  // assign target value as initial guess for injectors
1064  // only handles single phase control at the moment
1065  assert(numPhasesWithTargetsUnderThisControl == 1);
1066 
1067  for (int phase = 0; phase < np; ++phase) {
1068  if (distr[phase] > 0.) {
1069  xw.wellRates()[np*well_index + phase] = target / distr[phase];
1070  } else {
1071  xw.wellRates()[np * well_index + phase] = 0.;
1072  }
1073  }
1074  } else if (well_type_ == PRODUCER) {
1075  // update the rates of phases under control based on the target,
1076  // and also update rates of phases not under control to keep the rate ratio,
1077  // assuming the mobility ratio does not change for the production wells
1078  double original_rates_under_phase_control = 0.0;
1079  for (int phase = 0; phase < np; ++phase) {
1080  if (distr[phase] > 0.0) {
1081  original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase];
1082  }
1083  }
1084 
1085  if (original_rates_under_phase_control != 0.0 ) {
1086  double scaling_factor = target / original_rates_under_phase_control;
1087 
1088  for (int phase = 0; phase < np; ++phase) {
1089  xw.wellRates()[np * well_index + phase] *= scaling_factor;
1090  }
1091  } else { // scaling factor is not well defied when original_rates_under_phase_control is zero
1092  // separating targets equally between phases under control
1093  const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
1094  for (int phase = 0; phase < np; ++phase) {
1095  if (distr[phase] > 0.0) {
1096  xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase];
1097  } else {
1098  // this only happens for SURFACE_RATE control
1099  xw.wellRates()[np * well_index + phase] = target_rate_divided;
1100  }
1101  }
1102  }
1103  } else {
1104  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
1105  }
1106 
1107  break;
1108  } // end of switch
1109 
1110  updatePrimaryVariables(xw);
1111  }
1112 
1113 
1114 
1115 
1116 
1117  template<typename TypeTag>
1118  void
1120  computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
1121  const WellState& xw,
1122  std::vector<double>& b_perf,
1123  std::vector<double>& rsmax_perf,
1124  std::vector<double>& rvmax_perf,
1125  std::vector<double>& surf_dens_perf) const
1126  {
1127  const int nperf = number_of_perforations_;
1128  // TODO: can make this a member?
1129  const int numComp = numComponents();
1130  const PhaseUsage& pu = phaseUsage();
1131  b_perf.resize(nperf*numComp);
1132  surf_dens_perf.resize(nperf*numComp);
1133  const int w = index_of_well_;
1134 
1135  //rs and rv are only used if both oil and gas is present
1136  if (pu.phase_used[BlackoilPhases::Vapour] && pu.phase_used[BlackoilPhases::Liquid]) {
1137  rsmax_perf.resize(nperf);
1138  rvmax_perf.resize(nperf);
1139  }
1140 
1141  // Compute the average pressure in each well block
1142  for (int perf = 0; perf < nperf; ++perf) {
1143  const int cell_idx = well_cells_[perf];
1144  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
1145  const auto& fs = intQuants.fluidState();
1146 
1147  // TODO: this is another place to show why WellState need to be a vector of WellState.
1148  // TODO: to check why should be perf - 1
1149  const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1];
1150  const double p_avg = (xw.perfPress()[first_perf_ + perf] + p_above)/2;
1151  const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
1152 
1153  if (pu.phase_used[BlackoilPhases::Aqua]) {
1154  b_perf[ pu.phase_pos[BlackoilPhases::Aqua] + perf * numComp] =
1155  FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1156  }
1157 
1158  if (pu.phase_used[BlackoilPhases::Vapour]) {
1159  const int gaspos = pu.phase_pos[BlackoilPhases::Vapour] + perf * numComp;
1160  const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
1161 
1162  if (pu.phase_used[BlackoilPhases::Liquid]) {
1163  const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
1164  const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers
1165  rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
1166  if (oilrate > 0) {
1167  const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
1168  double rv = 0.0;
1169  if (gasrate > 0) {
1170  rv = oilrate / gasrate;
1171  }
1172  rv = std::min(rv, rvmax_perf[perf]);
1173 
1174  b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv);
1175  }
1176  else {
1177  b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1178  }
1179 
1180  } else {
1181  b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1182  }
1183  }
1184 
1185  if (pu.phase_used[BlackoilPhases::Liquid]) {
1186  const int oilpos = pu.phase_pos[BlackoilPhases::Liquid] + perf * numComp;
1187  const int oilpos_well = pu.phase_pos[BlackoilPhases::Liquid] + w * pu.num_phases;
1188  if (pu.phase_used[BlackoilPhases::Vapour]) {
1189  rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
1190  const int gaspos_well = pu.phase_pos[BlackoilPhases::Vapour] + w * pu.num_phases;
1191  const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
1192  if (gasrate > 0) {
1193  const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
1194  double rs = 0.0;
1195  if (oilrate > 0) {
1196  rs = gasrate / oilrate;
1197  }
1198  rs = std::min(rs, rsmax_perf[perf]);
1199  b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs);
1200  } else {
1201  b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1202  }
1203  } else {
1204  b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg);
1205  }
1206  }
1207 
1208  // Surface density.
1209  for (int p = 0; p < pu.num_phases; ++p) {
1210  surf_dens_perf[numComp*perf + p] = FluidSystem::referenceDensity( flowPhaseToEbosPhaseIdx( p ), fs.pvtRegionIndex());
1211  }
1212 
1213  // We use cell values for solvent injector
1214  if (has_solvent) {
1215  b_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value();
1216  surf_dens_perf[numComp*perf + contiSolventEqIdx] = intQuants.solventRefDensity();
1217  }
1218  }
1219  }
1220 
1221 
1222 
1223 
1224 
1225  template<typename TypeTag>
1226  void
1227  StandardWell<TypeTag>::
1228  computeConnectionDensities(const std::vector<double>& perfComponentRates,
1229  const std::vector<double>& b_perf,
1230  const std::vector<double>& rsmax_perf,
1231  const std::vector<double>& rvmax_perf,
1232  const std::vector<double>& surf_dens_perf)
1233  {
1234  // Verify that we have consistent input.
1235  const int np = number_of_phases_;
1236  const int nperf = number_of_perforations_;
1237  const int num_comp = numComponents();
1238  const PhaseUsage& phase_usage = phaseUsage();
1239 
1240  // 1. Compute the flow (in surface volume units for each
1241  // component) exiting up the wellbore from each perforation,
1242  // taking into account flow from lower in the well, and
1243  // in/out-flow at each perforation.
1244  std::vector<double> q_out_perf(nperf*num_comp);
1245 
1246  // TODO: investigate whether we should use the following techniques to calcuate the composition of flows in the wellbore
1247  // Iterate over well perforations from bottom to top.
1248  for (int perf = nperf - 1; perf >= 0; --perf) {
1249  for (int component = 0; component < num_comp; ++component) {
1250  if (perf == nperf - 1) {
1251  // This is the bottom perforation. No flow from below.
1252  q_out_perf[perf*num_comp+ component] = 0.0;
1253  } else {
1254  // Set equal to flow from below.
1255  q_out_perf[perf*num_comp + component] = q_out_perf[(perf+1)*num_comp + component];
1256  }
1257  // Subtract outflow through perforation.
1258  q_out_perf[perf*num_comp + component] -= perfComponentRates[perf*num_comp + component];
1259  }
1260  }
1261 
1262  // 2. Compute the component mix at each perforation as the
1263  // absolute values of the surface rates divided by their sum.
1264  // Then compute volume ratios (formation factors) for each perforation.
1265  // Finally compute densities for the segments associated with each perforation.
1266  const int gaspos = phase_usage.phase_pos[BlackoilPhases::Vapour];
1267  const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid];
1268  std::vector<double> mix(num_comp,0.0);
1269  std::vector<double> x(num_comp);
1270  std::vector<double> surf_dens(num_comp);
1271  std::vector<double> dens(nperf);
1272 
1273  for (int perf = 0; perf < nperf; ++perf) {
1274  // Find component mix.
1275  const double tot_surf_rate = std::accumulate(q_out_perf.begin() + num_comp*perf,
1276  q_out_perf.begin() + num_comp*(perf+1), 0.0);
1277  if (tot_surf_rate != 0.0) {
1278  for (int component = 0; component < num_comp; ++component) {
1279  mix[component] = std::fabs(q_out_perf[perf*num_comp + component]/tot_surf_rate);
1280  }
1281  } else {
1282  // No flow => use well specified fractions for mix.
1283  for (int phase = 0; phase < np; ++phase) {
1284  mix[phase] = comp_frac_[phase];
1285  }
1286  // intialize 0.0 for comIdx >= np;
1287  }
1288  // Compute volume ratio.
1289  x = mix;
1290  double rs = 0.0;
1291  double rv = 0.0;
1292  if (!rsmax_perf.empty() && mix[oilpos] > 0.0) {
1293  rs = std::min(mix[gaspos]/mix[oilpos], rsmax_perf[perf]);
1294  }
1295  if (!rvmax_perf.empty() && mix[gaspos] > 0.0) {
1296  rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
1297  }
1298  if (rs != 0.0) {
1299  // Subtract gas in oil from gas mixture
1300  x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
1301  }
1302  if (rv != 0.0) {
1303  // Subtract oil in gas from oil mixture
1304  x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
1305  }
1306  double volrat = 0.0;
1307  for (int component = 0; component < num_comp; ++component) {
1308  volrat += x[component] / b_perf[perf*num_comp+ component];
1309  }
1310  for (int component = 0; component < num_comp; ++component) {
1311  surf_dens[component] = surf_dens_perf[perf*num_comp+ component];
1312  }
1313 
1314  // Compute segment density.
1315  perf_densities_[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat;
1316  }
1317  }
1318 
1319 
1320 
1321 
1322  template<typename TypeTag>
1323  void
1324  StandardWell<TypeTag>::
1325  computeConnectionPressureDelta()
1326  {
1327  // Algorithm:
1328 
1329  // We'll assume the perforations are given in order from top to
1330  // bottom for each well. By top and bottom we do not necessarily
1331  // mean in a geometric sense (depth), but in a topological sense:
1332  // the 'top' perforation is nearest to the surface topologically.
1333  // Our goal is to compute a pressure delta for each perforation.
1334 
1335  // 1. Compute pressure differences between perforations.
1336  // dp_perf will contain the pressure difference between a
1337  // perforation and the one above it, except for the first
1338  // perforation for each well, for which it will be the
1339  // difference to the reference (bhp) depth.
1340 
1341  const int nperf = number_of_perforations_;
1342  perf_pressure_diffs_.resize(nperf, 0.0);
1343 
1344  for (int perf = 0; perf < nperf; ++perf) {
1345  const double z_above = perf == 0 ? ref_depth_ : perf_depth_[perf - 1];
1346  const double dz = perf_depth_[perf] - z_above;
1347  perf_pressure_diffs_[perf] = dz * perf_densities_[perf] * gravity_;
1348  }
1349 
1350  // 2. Compute pressure differences to the reference point (bhp) by
1351  // accumulating the already computed adjacent pressure
1352  // differences, storing the result in dp_perf.
1353  // This accumulation must be done per well.
1354  const auto beg = perf_pressure_diffs_.begin();
1355  const auto end = perf_pressure_diffs_.end();
1356  std::partial_sum(beg, end, beg);
1357  }
1358 
1359 
1360 
1361 
1362 
1363  template<typename TypeTag>
1364  typename StandardWell<TypeTag>::ConvergenceReport
1366  getWellConvergence(const std::vector<double>& B_avg) const
1367  {
1368  typedef double Scalar;
1369  typedef std::vector< Scalar > Vector;
1370 
1371  const int np = number_of_phases_;
1372  const int numComp = numComponents();
1373 
1374  // the following implementation assume that the polymer is always after the w-o-g phases
1375  // For the polymer case, there is one more mass balance equations of reservoir than wells
1376  assert((int(B_avg.size()) == numComp) || has_polymer);
1377 
1378  const double tol_wells = param_.tolerance_wells_;
1379  const double maxResidualAllowed = param_.max_residual_allowed_;
1380 
1381  // TODO: it should be the number of numWellEq
1382  // using numComp here for flow_ebos running 2p case.
1383  std::vector<Scalar> res(numComp);
1384  for (int comp = 0; comp < numComp; ++comp) {
1385  // magnitude of the residual matters
1386  res[comp] = std::abs(resWell_[0][comp]);
1387  }
1388 
1389  Vector well_flux_residual(numComp);
1390 
1391  // Finish computation
1392  for ( int compIdx = 0; compIdx < numComp; ++compIdx )
1393  {
1394  well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx];
1395  }
1396 
1397  ConvergenceReport report;
1398  // checking if any NaN or too large residuals found
1399  // TODO: not understand why phase here while component in other places.
1400  for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
1401  const auto& phaseName = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(phaseIdx));
1402 
1403  if (std::isnan(well_flux_residual[phaseIdx])) {
1404  report.nan_residual_found = true;
1405  const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName};
1406  report.nan_residual_wells.push_back(problem_well);
1407  } else {
1408  if (well_flux_residual[phaseIdx] > maxResidualAllowed) {
1409  report.too_large_residual_found = true;
1410  const typename ConvergenceReport::ProblemWell problem_well = {name(), phaseName};
1411  report.too_large_residual_wells.push_back(problem_well);
1412  }
1413  }
1414  }
1415 
1416  if ( !(report.nan_residual_found || report.too_large_residual_found) ) { // no abnormal residual value found
1417  // check convergence
1418  for ( int compIdx = 0; compIdx < numComp; ++compIdx )
1419  {
1420  report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells);
1421  }
1422  } else { // abnormal values found and no need to check the convergence
1423  report.converged = false;
1424  }
1425 
1426  return report;
1427  }
1428 
1429 
1430 
1431 
1432 
1433  template<typename TypeTag>
1434  void
1436  computeWellConnectionDensitesPressures(const WellState& xw,
1437  const std::vector<double>& b_perf,
1438  const std::vector<double>& rsmax_perf,
1439  const std::vector<double>& rvmax_perf,
1440  const std::vector<double>& surf_dens_perf)
1441  {
1442  // Compute densities
1443  const int nperf = number_of_perforations_;
1444  const int numComponent = numComponents();
1445  const int np = number_of_phases_;
1446  std::vector<double> perfRates(b_perf.size(),0.0);
1447 
1448  for (int perf = 0; perf < nperf; ++perf) {
1449  for (int phase = 0; phase < np; ++phase) {
1450  perfRates[perf*numComponent + phase] = xw.perfPhaseRates()[(first_perf_ + perf) * np + phase];
1451  }
1452  if(has_solvent) {
1453  perfRates[perf*numComponent + contiSolventEqIdx] = xw.perfRateSolvent()[first_perf_ + perf];
1454  }
1455  }
1456 
1457  computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1458 
1459  computeConnectionPressureDelta();
1460 
1461  }
1462 
1463 
1464 
1465 
1466 
1467  template<typename TypeTag>
1468  void
1469  StandardWell<TypeTag>::
1470  computeWellConnectionPressures(const Simulator& ebosSimulator,
1471  const WellState& well_state)
1472  {
1473  // 1. Compute properties required by computeConnectionPressureDelta().
1474  // Note that some of the complexity of this part is due to the function
1475  // taking std::vector<double> arguments, and not Eigen objects.
1476  std::vector<double> b_perf;
1477  std::vector<double> rsmax_perf;
1478  std::vector<double> rvmax_perf;
1479  std::vector<double> surf_dens_perf;
1480  computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1481  computeWellConnectionDensitesPressures(well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1482  }
1483 
1484 
1485 
1486 
1487 
1488  template<typename TypeTag>
1489  void
1490  StandardWell<TypeTag>::
1491  solveEqAndUpdateWellState(WellState& well_state)
1492  {
1493  // We assemble the well equations, then we check the convergence,
1494  // which is why we do not put the assembleWellEq here.
1495  BVectorWell dx_well(1);
1496  invDuneD_.mv(resWell_, dx_well);
1497 
1498  updateWellState(dx_well, well_state);
1499  }
1500 
1501 
1502 
1503 
1504 
1505  template<typename TypeTag>
1506  void
1507  StandardWell<TypeTag>::
1508  calculateExplicitQuantities(const Simulator& ebosSimulator,
1509  const WellState& well_state)
1510  {
1511  computeWellConnectionPressures(ebosSimulator, well_state);
1512  computeAccumWell();
1513  }
1514 
1515 
1516 
1517 
1518 
1519  template<typename TypeTag>
1520  void
1521  StandardWell<TypeTag>::
1522  computeAccumWell()
1523  {
1524  // TODO: it should be num_comp, while it also bring problem for
1525  // the polymer case.
1526  for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
1527  F0_[eq_idx] = wellSurfaceVolumeFraction(eq_idx).value();
1528  }
1529  }
1530 
1531 
1532 
1533 
1534 
1535  template<typename TypeTag>
1536  void
1538  apply(const BVector& x, BVector& Ax) const
1539  {
1540  assert( Bx_.size() == duneB_.N() );
1541  assert( invDrw_.size() == invDuneD_.N() );
1542 
1543  // Bx_ = duneB_ * x
1544  duneB_.mv(x, Bx_);
1545  // invDBx = invDuneD_ * Bx_
1546  // TODO: with this, we modified the content of the invDrw_.
1547  // Is it necessary to do this to save some memory?
1548  BVectorWell& invDBx = invDrw_;
1549  invDuneD_.mv(Bx_, invDBx);
1550 
1551  // Ax = Ax - duneC_^T * invDBx
1552  duneC_.mmtv(invDBx,Ax);
1553  }
1554 
1555 
1556 
1557 
1558  template<typename TypeTag>
1559  void
1561  apply(BVector& r) const
1562  {
1563  assert( invDrw_.size() == invDuneD_.N() );
1564 
1565  // invDrw_ = invDuneD_ * resWell_
1566  invDuneD_.mv(resWell_, invDrw_);
1567  // r = r - duneC_^T * invDrw_
1568  duneC_.mmtv(invDrw_, r);
1569  }
1570 
1571 
1572 
1573 
1574 
1575  template<typename TypeTag>
1576  void
1578  recoverSolutionWell(const BVector& x, BVectorWell& xw) const
1579  {
1580  BVectorWell resWell = resWell_;
1581  // resWell = resWell - B * x
1582  duneB_.mmv(x, resWell);
1583  // xw = D^-1 * resWell
1584  invDuneD_.mv(resWell, xw);
1585  }
1586 
1587 
1588 
1589 
1590 
1591  template<typename TypeTag>
1592  void
1595  WellState& well_state) const
1596  {
1597  BVectorWell xw(1);
1598  recoverSolutionWell(x, xw);
1599  updateWellState(xw, well_state);
1600  }
1601 
1602 
1603 
1604 
1605 
1606  template<typename TypeTag>
1607  void
1609  computeWellRatesWithBhp(const Simulator& ebosSimulator,
1610  const EvalWell& bhp,
1611  std::vector<double>& well_flux) const
1612  {
1613  const int np = number_of_phases_;
1614  const int numComp = numComponents();
1615  well_flux.resize(np, 0.0);
1616 
1617  const bool allow_cf = crossFlowAllowed(ebosSimulator);
1618 
1619  for (int perf = 0; perf < number_of_perforations_; ++perf) {
1620  const int cell_idx = well_cells_[perf];
1621  const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
1622  // flux for each perforation
1623  std::vector<EvalWell> cq_s(numComp, 0.0);
1624  std::vector<EvalWell> mob(numComp, 0.0);
1625  getMobility(ebosSimulator, perf, mob);
1626  computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s);
1627 
1628  for(int p = 0; p < np; ++p) {
1629  well_flux[p] += cq_s[p].value();
1630  }
1631  }
1632  }
1633 
1634 
1635 
1636 
1637 
1638  template<typename TypeTag>
1639  std::vector<double>
1640  StandardWell<TypeTag>::
1641  computeWellPotentialWithTHP(const Simulator& ebosSimulator,
1642  const double initial_bhp, // bhp from BHP constraints
1643  const std::vector<double>& initial_potential) const
1644  {
1645  // TODO: pay attention to the situation that finally the potential is calculated based on the bhp control
1646  // TODO: should we consider the bhp constraints during the iterative process?
1647  const int np = number_of_phases_;
1648 
1649  assert( np == int(initial_potential.size()) );
1650 
1651  std::vector<double> potentials = initial_potential;
1652  std::vector<double> old_potentials = potentials; // keeping track of the old potentials
1653 
1654  double bhp = initial_bhp;
1655  double old_bhp = bhp;
1656 
1657  bool converged = false;
1658  const int max_iteration = 1000;
1659  const double bhp_tolerance = 1000.; // 1000 pascal
1660 
1661  int iteration = 0;
1662 
1663  while ( !converged && iteration < max_iteration ) {
1664  // for each iteration, we calculate the bhp based on the rates/potentials with thp constraints
1665  // with considering the bhp value from the bhp limits. At the beginning of each iteration,
1666  // we initialize the bhp to be the bhp value from the bhp limits. Then based on the bhp values calculated
1667  // from the thp constraints, we decide the effective bhp value for well potential calculation.
1668  bhp = initial_bhp;
1669 
1670  // The number of the well controls/constraints
1671  const int nwc = well_controls_get_num(well_controls_);
1672 
1673  for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
1674  if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
1675  const Opm::PhaseUsage& pu = phaseUsage();
1676 
1677  std::vector<double> rates(3, 0.0);
1678  if (active()[ Water ]) {
1679  rates[ Water ] = potentials[pu.phase_pos[ Water ] ];
1680  }
1681  if (active()[ Oil ]) {
1682  rates[ Oil ] = potentials[pu.phase_pos[ Oil ] ];
1683  }
1684  if (active()[ Gas ]) {
1685  rates[ Gas ] = potentials[pu.phase_pos[ Gas ] ];
1686  }
1687 
1688  const double bhp_calculated = calculateBhpFromThp(rates, ctrl_index);
1689 
1690  if (well_type_ == INJECTOR && bhp_calculated < bhp ) {
1691  bhp = bhp_calculated;
1692  }
1693 
1694  if (well_type_ == PRODUCER && bhp_calculated > bhp) {
1695  bhp = bhp_calculated;
1696  }
1697  }
1698  }
1699 
1700  // there should be always some available bhp/thp constraints there
1701  if (std::isinf(bhp) || std::isnan(bhp)) {
1702  OPM_THROW(std::runtime_error, "Unvalid bhp value obtained during the potential calculation for well " << name());
1703  }
1704 
1705  converged = std::abs(old_bhp - bhp) < bhp_tolerance;
1706 
1707  computeWellRatesWithBhp(ebosSimulator, bhp, potentials);
1708 
1709  // checking whether the potentials have valid values
1710  for (const double value : potentials) {
1711  if (std::isinf(value) || std::isnan(value)) {
1712  OPM_THROW(std::runtime_error, "Unvalid potential value obtained during the potential calculation for well " << name());
1713  }
1714  }
1715 
1716  if (!converged) {
1717  old_bhp = bhp;
1718  for (int p = 0; p < np; ++p) {
1719  // TODO: improve the interpolation, will it always be valid with the way below?
1720  // TODO: finding better paramters, better iteration strategy for better convergence rate.
1721  const double potential_update_damping_factor = 0.001;
1722  potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p];
1723  old_potentials[p] = potentials[p];
1724  }
1725  }
1726 
1727  ++iteration;
1728  }
1729 
1730  if (!converged) {
1731  OPM_THROW(std::runtime_error, "Failed in getting converged for the potential calculation for well " << name());
1732  }
1733 
1734  return potentials;
1735  }
1736 
1737 
1738 
1739 
1740 
1741  template<typename TypeTag>
1742  void
1744  computeWellPotentials(const Simulator& ebosSimulator,
1745  const WellState& well_state,
1746  std::vector<double>& well_potentials) // const
1747  {
1748  updatePrimaryVariables(well_state);
1749  computeWellConnectionPressures(ebosSimulator, well_state);
1750 
1751  // initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
1752  // TODO: for computeWellPotentials, no derivative is required actually
1753  initPrimaryVariablesEvaluation();
1754 
1755  const int np = number_of_phases_;
1756  well_potentials.resize(np, 0.0);
1757 
1758  // get the bhp value based on the bhp constraints
1759  const double bhp = mostStrictBhpFromBhpLimits();
1760 
1761  // does the well have a THP related constraint?
1762  if ( !wellHasTHPConstraints() ) {
1763  assert(std::abs(bhp) != std::numeric_limits<double>::max());
1764 
1765  computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
1766  } else {
1767  // the well has a THP related constraint
1768  // checking whether a well is newly added, it only happens at the beginning of the report step
1769  if ( !well_state.isNewWell(index_of_well_) ) {
1770  for (int p = 0; p < np; ++p) {
1771  // This is dangerous for new added well
1772  // since we are not handling the initialization correctly for now
1773  well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p];
1774  }
1775  } else {
1776  // We need to generate a reasonable rates to start the iteration process
1777  computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
1778  for (double& value : well_potentials) {
1779  // make the value a little safer in case the BHP limits are default ones
1780  // TODO: a better way should be a better rescaling based on the investigation of the VFP table.
1781  const double rate_safety_scaling_factor = 0.00001;
1782  value *= rate_safety_scaling_factor;
1783  }
1784  }
1785 
1786  well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials);
1787  }
1788  }
1789 
1790 
1791 
1792 
1793 
1794  template<typename TypeTag>
1795  void
1797  updatePrimaryVariables(const WellState& well_state) const
1798  {
1799  const int np = number_of_phases_;
1800  const int well_index = index_of_well_;
1801  const WellControls* wc = well_controls_;
1802  const double* distr = well_controls_get_current_distr(wc);
1803  const auto pu = phaseUsage();
1804 
1805  switch (well_controls_get_current_type(wc)) {
1806  case THP:
1807  case BHP: {
1808  primary_variables_[XvarWell] = 0.0;
1809  if (well_type_ == INJECTOR) {
1810  for (int p = 0; p < np; ++p) {
1811  primary_variables_[XvarWell] += well_state.wellRates()[np*well_index + p] * comp_frac_[p];
1812  }
1813  } else {
1814  for (int p = 0; p < np; ++p) {
1815  primary_variables_[XvarWell] += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
1816  }
1817  }
1818  break;
1819  }
1820  case RESERVOIR_RATE: // Intentional fall-through
1821  case SURFACE_RATE:
1822  primary_variables_[XvarWell] = well_state.bhp()[well_index];
1823  break;
1824  } // end of switch
1825 
1826  double tot_well_rate = 0.0;
1827  for (int p = 0; p < np; ++p) {
1828  tot_well_rate += scalingFactor(p) * well_state.wellRates()[np*well_index + p];
1829  }
1830  if(std::abs(tot_well_rate) > 0) {
1831  if (active()[ Water ]) {
1832  primary_variables_[WFrac] = scalingFactor(pu.phase_pos[Water]) * well_state.wellRates()[np*well_index + pu.phase_pos[Water]] / tot_well_rate;
1833  }
1834  if (active()[ Gas ]) {
1835  primary_variables_[GFrac] = scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates()[np*well_index + pu.phase_pos[Gas]] - well_state.solventWellRate(well_index)) / tot_well_rate ;
1836  }
1837  if (has_solvent) {
1838  primary_variables_[SFrac] = scalingFactor(pu.phase_pos[Gas]) * well_state.solventWellRate(well_index) / tot_well_rate ;
1839  }
1840  } else { // tot_well_rate == 0
1841  if (well_type_ == INJECTOR) {
1842  // only single phase injection handled
1843  if (active()[Water]) {
1844  if (distr[Water] > 0.0) {
1845  primary_variables_[WFrac] = 1.0;
1846  } else {
1847  primary_variables_[WFrac] = 0.0;
1848  }
1849  }
1850 
1851  if (active()[Gas]) {
1852  if (distr[pu.phase_pos[Gas]] > 0.0) {
1853  primary_variables_[GFrac] = 1.0 - wsolvent();
1854  if (has_solvent) {
1855  primary_variables_[SFrac] = wsolvent();
1856  }
1857  } else {
1858  primary_variables_[GFrac] = 0.0;
1859  }
1860  }
1861 
1862  // TODO: it is possible to leave injector as a oil well,
1863  // when F_w and F_g both equals to zero, not sure under what kind of circumstance
1864  // this will happen.
1865  } else if (well_type_ == PRODUCER) { // producers
1866  // TODO: the following are not addressed for the solvent case yet
1867  if (active()[Water]) {
1868  primary_variables_[WFrac] = 1.0 / np;
1869  }
1870  if (active()[Gas]) {
1871  primary_variables_[GFrac] = 1.0 / np;
1872  }
1873  } else {
1874  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
1875  }
1876  }
1877  }
1878 
1879 
1880 
1881 
1882 
1883  template<typename TypeTag>
1884  template<class ValueType>
1885  ValueType
1886  StandardWell<TypeTag>::
1887  calculateBhpFromThp(const std::vector<ValueType>& rates,
1888  const int control_index) const
1889  {
1890  // TODO: when well is under THP control, the BHP is dependent on the rates,
1891  // the well rates is also dependent on the BHP, so it might need to do some iteration.
1892  // However, when group control is involved, change of the rates might impacts other wells
1893  // so iterations on a higher level will be required. Some investigation might be needed when
1894  // we face problems under THP control.
1895 
1896  assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
1897 
1898  const ValueType aqua = rates[Water];
1899  const ValueType liquid = rates[Oil];
1900  const ValueType vapour = rates[Gas];
1901 
1902  const int vfp = well_controls_iget_vfp(well_controls_, control_index);
1903  const double& thp = well_controls_iget_target(well_controls_, control_index);
1904  const double& alq = well_controls_iget_alq(well_controls_, control_index);
1905 
1906  // pick the density in the top layer
1907  const double rho = perf_densities_[0];
1908 
1909  ValueType bhp = 0.;
1910  if (well_type_ == INJECTOR) {
1911  const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
1912 
1913  const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1914 
1915  bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
1916  }
1917  else if (well_type_ == PRODUCER) {
1918  const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
1919 
1920  const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1921 
1922  bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
1923  }
1924  else {
1925  OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
1926  }
1927 
1928  return bhp;
1929  }
1930 
1931 
1932 
1933 
1934 
1935  template<typename TypeTag>
1936  double
1937  StandardWell<TypeTag>::
1938  calculateThpFromBhp(const std::vector<double>& rates,
1939  const int control_index,
1940  const double bhp) const
1941  {
1942  assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
1943 
1944  const double aqua = rates[Water];
1945  const double liquid = rates[Oil];
1946  const double vapour = rates[Gas];
1947 
1948  const int vfp = well_controls_iget_vfp(well_controls_, control_index);
1949  const double& alq = well_controls_iget_alq(well_controls_, control_index);
1950 
1951  // pick the density in the top layer
1952  const double rho = perf_densities_[0];
1953 
1954  double thp = 0.0;
1955  if (well_type_ == INJECTOR) {
1956  const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
1957 
1958  const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1959 
1960  thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp);
1961  }
1962  else if (well_type_ == PRODUCER) {
1963  const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
1964 
1965  const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
1966 
1967  thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq);
1968  }
1969  else {
1970  OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
1971  }
1972 
1973  return thp;
1974  }
1975 
1976  template<typename TypeTag>
1977  double
1978  StandardWell<TypeTag>::scalingFactor(const int phaseIdx) const
1979  {
1980  const WellControls* wc = well_controls_;
1981  const double* distr = well_controls_get_current_distr(wc);
1982 
1983  if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
1984  if (has_solvent && phaseIdx == contiSolventEqIdx )
1985  OPM_THROW(std::runtime_error, "RESERVOIR_RATE control in combination with solvent is not implemented");
1986 
1987  return distr[phaseIdx];
1988  }
1989  const auto& pu = phaseUsage();
1990  if (active()[Water] && pu.phase_pos[Water] == phaseIdx)
1991  return 1.0;
1992  if (active()[Oil] && pu.phase_pos[Oil] == phaseIdx)
1993  return 1.0;
1994  if (active()[Gas] && pu.phase_pos[Gas] == phaseIdx)
1995  return 0.01;
1996  if (has_solvent && phaseIdx == contiSolventEqIdx )
1997  return 0.01;
1998 
1999  // we should not come this far
2000  assert(false);
2001  return 1.0;
2002  }
2003 
2004 
2005 }
virtual ConvergenceReport getWellConvergence(const std::vector< double > &B_avg) const
check whether the well equations get converged for this well
Definition: StandardWell_impl.hpp:1366
Definition: StandardWell.hpp:33
a struct to collect information about the convergence checking
Definition: WellInterface.hpp:117
virtual void apply(const BVector &x, BVector &Ax) const
Ax = Ax - C D^-1 B x.
Definition: StandardWell_impl.hpp:1538
virtual void updateWellStateWithTarget(const int current, WellState &xw) const
updating the well state based the control mode specified with current
Definition: StandardWell_impl.hpp:1013
virtual void computeWellPotentials(const Simulator &ebosSimulator, const WellState &well_state, std::vector< double > &well_potentials)
computing the well potentials for group control
Definition: StandardWell_impl.hpp:1744
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
virtual void recoverWellSolutionAndUpdateWellState(const BVector &x, WellState &well_state) const
using the solution x to recover the solution xw for wells and applying xw to update Well State ...
Definition: StandardWell_impl.hpp:1594