All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
StandardWells_impl.hpp
1 /*
2  Copyright 2016 SINTEF ICT, Applied Mathematics.
3  Copyright 2016 Statoil ASA.
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #include <opm/autodiff/StandardWells.hpp>
23 #include <opm/autodiff/WellDensitySegmented.hpp>
24 
25 #include <opm/autodiff/VFPInjProperties.hpp>
26 #include <opm/autodiff/VFPProdProperties.hpp>
27 #include <opm/autodiff/WellHelpers.hpp>
28 
29 
30 
31 
32 namespace Opm
33 {
34 
35 
36  StandardWells::
37  WellOps::WellOps(const Wells* wells)
38  : w2p(),
39  p2w(),
40  well_cells()
41  {
42  if( wells )
43  {
44  w2p = Eigen::SparseMatrix<double>(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
45  p2w = Eigen::SparseMatrix<double>(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
46 
47  const int nw = wells->number_of_wells;
48  const int* const wpos = wells->well_connpos;
49 
50  typedef Eigen::Triplet<double> Tri;
51 
52  std::vector<Tri> scatter, gather;
53  scatter.reserve(wpos[nw]);
54  gather .reserve(wpos[nw]);
55 
56  for (int w = 0, i = 0; w < nw; ++w) {
57  for (; i < wpos[ w + 1 ]; ++i) {
58  scatter.push_back(Tri(i, w, 1.0));
59  gather .push_back(Tri(w, i, 1.0));
60  }
61  }
62 
63  w2p.setFromTriplets(scatter.begin(), scatter.end());
64  p2w.setFromTriplets(gather .begin(), gather .end());
65 
66  well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]);
67  }
68  }
69 
70 
71 
72 
73 
74  StandardWells::StandardWells(const Wells* wells_arg, WellCollection* well_collection)
75  : wells_active_(wells_arg!=nullptr)
76  , wells_(wells_arg)
77  , wops_(wells_arg)
78  , well_collection_(well_collection)
79  , well_perforation_efficiency_factors_(Vector::Ones(wells_!=nullptr ? wells_->well_connpos[wells_->number_of_wells] : 0))
80  , fluid_(nullptr)
81  , active_(nullptr)
82  , phase_condition_(nullptr)
83  , vfp_properties_(nullptr)
84  , well_perforation_densities_(Vector())
85  , well_perforation_pressure_diffs_(Vector())
86  , store_well_perforation_fluxes_(false)
87  {
88  }
89 
90 
91 
92 
93 
94  void
95  StandardWells::init(const BlackoilPropsAdFromDeck* fluid_arg,
96  const std::vector<bool>* active_arg,
97  const std::vector<PhasePresence>* pc_arg,
98  const VFPProperties* vfp_properties_arg,
99  const double gravity_arg,
100  const Vector& depth_arg)
101  {
102  fluid_ = fluid_arg;
103  active_ = active_arg;
104  phase_condition_ = pc_arg;
105  vfp_properties_ = vfp_properties_arg;
106  gravity_ = gravity_arg;
107  perf_cell_depth_ = subset(depth_arg, wellOps().well_cells);;
108 
109  calculateEfficiencyFactors();
110  }
111 
112 
113 
114 
115 
116  const Wells& StandardWells::wells() const
117  {
118  assert(wells_ != 0);
119  return *(wells_);
120  }
121 
122 
123  const Wells* StandardWells::wellsPointer() const
124  {
125  return wells_;
126  }
127 
128 
129 
131  {
132  return wells_active_;
133  }
134 
135 
136 
137 
138 
139  void StandardWells::setWellsActive(const bool wells_active)
140  {
141  wells_active_ = wells_active;
142  }
143 
144 
145 
146 
147 
149  {
150  return wells_ ? (wells_->number_of_wells > 0 ) : false;
151  }
152 
153 
154 
155 
156 
157  int
158  StandardWells::numWellVars() const
159  {
160  if ( !localWellsActive() )
161  {
162  return 0;
163  }
164 
165  // For each well, we have a bhp variable, and one flux per phase.
166  const int nw = wells().number_of_wells;
167  return (numPhases() + 1) * nw;
168  }
169 
170 
171 
172 
173 
174  const StandardWells::WellOps&
175  StandardWells::wellOps() const
176  {
177  return wops_;
178  }
179 
180 
181 
182 
183 
185  {
186  return well_perforation_densities_;
187  }
188 
189 
190 
191 
192 
193  const StandardWells::Vector&
195  {
196  return well_perforation_densities_;
197  }
198 
199 
200 
201 
202 
203  StandardWells::Vector&
205  {
206  return well_perforation_pressure_diffs_;
207  }
208 
209 
210 
211 
212 
213  const StandardWells::Vector&
215  {
216  return well_perforation_pressure_diffs_;
217  }
218 
219 
220 
221 
222  template<class SolutionState, class WellState>
223  void
224  StandardWells::
225  computePropertiesForWellConnectionPressures(const SolutionState& state,
226  const WellState& xw,
227  std::vector<double>& b_perf,
228  std::vector<double>& rsmax_perf,
229  std::vector<double>& rvmax_perf,
230  std::vector<double>& surf_dens_perf)
231  {
232  const int nperf = wells().well_connpos[wells().number_of_wells];
233  const int nw = wells().number_of_wells;
234 
235  // Compute the average pressure in each well block
236  const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf);
237  Vector avg_press = perf_press*0;
238  for (int w = 0; w < nw; ++w) {
239  for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
240  const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
241  const double p_avg = (perf_press[perf] + p_above)/2;
242  avg_press[perf] = p_avg;
243  }
244  }
245 
246  const std::vector<int>& well_cells = wellOps().well_cells;
247 
248  // Use cell values for the temperature as the wells don't knows its temperature yet.
249  const ADB perf_temp = subset(state.temperature, well_cells);
250 
251  // Compute b, rsmax, rvmax values for perforations.
252  // Evaluate the properties using average well block pressures
253  // and cell values for rs, rv, phase condition and temperature.
254  const ADB avg_press_ad = ADB::constant(avg_press);
255  std::vector<PhasePresence> perf_cond(nperf);
256  // const std::vector<PhasePresence>& pc = phaseCondition();
257  for (int perf = 0; perf < nperf; ++perf) {
258  perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
259  }
260  const PhaseUsage& pu = fluid_->phaseUsage();
261  DataBlock b(nperf, pu.num_phases);
262  if (pu.phase_used[BlackoilPhases::Aqua]) {
263  const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
264  b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
265  }
266  assert((*active_)[Oil]);
267  const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
268  if (pu.phase_used[BlackoilPhases::Liquid]) {
269  const ADB perf_rs = (state.rs.size() > 0) ? subset(state.rs, well_cells) : ADB::null();
270  const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
271  b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
272  }
273  if (pu.phase_used[BlackoilPhases::Vapour]) {
274  const ADB perf_rv = (state.rv.size() > 0) ? subset(state.rv, well_cells) : ADB::null();
275  const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
276  b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
277  }
278  if (pu.phase_used[BlackoilPhases::Liquid] && pu.phase_used[BlackoilPhases::Vapour]) {
279  const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
280  rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
281 
282  const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
283  rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
284  }
285 
286  // b is row major, so can just copy data.
287  b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
288 
289  // Surface density.
290  // The compute density segment wants the surface densities as
291  // an np * number of wells cells array
292  Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
293  for (int phase = 1; phase < pu.num_phases; ++phase) {
294  rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
295  }
296  surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases);
297 
298  }
299 
300 
301 
302 
303 
304  template <class WellState>
305  void
306  StandardWells::
307  computeWellConnectionDensitesPressures(const WellState& xw,
308  const std::vector<double>& b_perf,
309  const std::vector<double>& rsmax_perf,
310  const std::vector<double>& rvmax_perf,
311  const std::vector<double>& surf_dens_perf,
312  const std::vector<double>& depth_perf,
313  const double grav)
314  {
315  // Compute densities
316  std::vector<double> cd =
318  wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
319  b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
320 
321  const int nperf = wells().well_connpos[wells().number_of_wells];
322 
323  // Compute pressure deltas
324  std::vector<double> cdp =
326  wells(), depth_perf, cd, grav);
327 
328  // Store the results
329  well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf);
330  well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf);
331  }
332 
333 
334 
335 
336 
337  template <class SolutionState, class WellState>
338  void
339  StandardWells::
340  computeWellConnectionPressures(const SolutionState& state,
341  const WellState& xw)
342  {
343  if( ! localWellsActive() ) return ;
344  // 1. Compute properties required by computeConnectionPressureDelta().
345  // Note that some of the complexity of this part is due to the function
346  // taking std::vector<double> arguments, and not Eigen objects.
347  std::vector<double> b_perf;
348  std::vector<double> rsmax_perf;
349  std::vector<double> rvmax_perf;
350  std::vector<double> surf_dens_perf;
351  computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
352 
353  const Vector& pdepth = perf_cell_depth_;
354  const int nperf = wells().well_connpos[wells().number_of_wells];
355  const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
356 
357  computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
358 
359  }
360 
361 
362 
363 
364 
365  template <class ReservoirResidualQuant, class SolutionState>
366  void
367  StandardWells::
368  extractWellPerfProperties(const SolutionState& /* state */,
369  const std::vector<ReservoirResidualQuant>& rq,
370  std::vector<ADB>& mob_perfcells,
371  std::vector<ADB>& b_perfcells) const
372  {
373  // If we have wells, extract the mobilities and b-factors for
374  // the well-perforated cells.
375  if ( !localWellsActive() ) {
376  mob_perfcells.clear();
377  b_perfcells.clear();
378  return;
379  } else {
380  const std::vector<int>& well_cells = wellOps().well_cells;
381  const int num_phases = wells().number_of_phases;
382  mob_perfcells.resize(num_phases, ADB::null());
383  b_perfcells.resize(num_phases, ADB::null());
384  for (int phase = 0; phase < num_phases; ++phase) {
385  mob_perfcells[phase] = subset(rq[phase].mob, well_cells);
386  b_perfcells[phase] = subset(rq[phase].b, well_cells);
387  }
388  }
389  }
390 
391 
392 
393 
394 
395 
396  template <class SolutionState>
397  void
398  StandardWells::
399  computeWellFlux(const SolutionState& state,
400  const std::vector<ADB>& mob_perfcells,
401  const std::vector<ADB>& b_perfcells,
402  Vector& aliveWells,
403  std::vector<ADB>& cq_s) const
404  {
405  if( ! localWellsActive() ) return ;
406 
407  const int np = wells().number_of_phases;
408  const int nw = wells().number_of_wells;
409  const int nperf = wells().well_connpos[nw];
410  Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
411  const std::vector<int>& well_cells = wellOps().well_cells;
412 
413  // pressure diffs computed already (once per step, not changing per iteration)
414  const Vector& cdp = wellPerforationPressureDiffs();
415  // Extract needed quantities for the perforation cells
416  const ADB& p_perfcells = subset(state.pressure, well_cells);
417 
418  // Perforation pressure
419  const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
420 
421  // Pressure drawdown (also used to determine direction of flow)
422  const ADB drawdown = p_perfcells - perfpressure;
423 
424  // Compute vectors with zero and ones that
425  // selects the wanted quantities.
426 
427  // selects injection perforations
428  Vector selectInjectingPerforations = Vector::Zero(nperf);
429  // selects producing perforations
430  Vector selectProducingPerforations = Vector::Zero(nperf);
431  for (int c = 0; c < nperf; ++c){
432  if (drawdown.value()[c] < 0)
433  selectInjectingPerforations[c] = 1;
434  else
435  selectProducingPerforations[c] = 1;
436  }
437 
438  // Handle cross flow
439  const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
440  const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
441  for (int w = 0; w < nw; ++w) {
442  if (!wells().allow_cf[w]) {
443  for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
444  // Crossflow is not allowed; reverse flow is prevented.
445  // At least one of the perforation must be open in order to have a meeningful
446  // equation to solve. For the special case where all perforations have reverse flow,
447  // and the target rate is non-zero all of the perforations are keept open.
448  if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
449  selectProducingPerforations[perf] = 0.0;
450  } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
451  selectInjectingPerforations[perf] = 0.0;
452  }
453  }
454  }
455  }
456 
457  // HANDLE FLOW INTO WELLBORE
458  // compute phase volumetric rates at standard conditions
459  std::vector<ADB> cq_p(np, ADB::null());
460  std::vector<ADB> cq_ps(np, ADB::null());
461  for (int phase = 0; phase < np; ++phase) {
462  cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
463  cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
464  }
465  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
466  if ((*active_)[Oil] && (*active_)[Gas]) {
467  const int oilpos = pu.phase_pos[Oil];
468  const int gaspos = pu.phase_pos[Gas];
469  const ADB cq_psOil = cq_ps[oilpos];
470  const ADB cq_psGas = cq_ps[gaspos];
471  const ADB& rv_perfcells = subset(state.rv, well_cells);
472  const ADB& rs_perfcells = subset(state.rs, well_cells);
473  cq_ps[gaspos] += rs_perfcells * cq_psOil;
474  cq_ps[oilpos] += rv_perfcells * cq_psGas;
475  }
476 
477  // HANDLE FLOW OUT FROM WELLBORE
478  // Using total mobilities
479  ADB total_mob = mob_perfcells[0];
480  for (int phase = 1; phase < np; ++phase) {
481  total_mob += mob_perfcells[phase];
482  }
483  // injection perforations total volume rates
484  const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
485 
486  // Store well perforation total fluxes (reservor volumes) if requested.
487  if (store_well_perforation_fluxes_) {
488  // Ugly const-cast, but unappealing alternatives.
489  Vector& wf = const_cast<Vector&>(well_perforation_fluxes_);
490  wf = cqt_i.value();
491  for (int phase = 0; phase < np; ++phase) {
492  wf += cq_p[phase].value();
493  }
494  }
495 
496  // compute wellbore mixture for injecting perforations
497  // The wellbore mixture depends on the inflow from the reservoar
498  // and the well injection rates.
499 
500  // compute avg. and total wellbore phase volumetric rates at standard conds
501  const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
502  std::vector<ADB> wbq(np, ADB::null());
503  ADB wbqt = ADB::constant(Vector::Zero(nw));
504  for (int phase = 0; phase < np; ++phase) {
505  const ADB& q_ps = wellOps().p2w * cq_ps[phase];
506  const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
507  Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
508  const int pos = pu.phase_pos[phase];
509  wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps;
510  wbqt += wbq[phase];
511  }
512  // compute wellbore mixture at standard conditions.
513  Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
514  std::vector<ADB> cmix_s(np, ADB::null());
515  for (int phase = 0; phase < np; ++phase) {
516  const int pos = pu.phase_pos[phase];
517  cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
518  }
519 
520  // compute volume ratio between connection at standard conditions
521  ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
522 
523  if ((*active_)[Water]) {
524  const int watpos = pu.phase_pos[Water];
525  volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
526  }
527 
528  if ((*active_)[Oil] && (*active_)[Gas]) {
529  // Incorporate RS/RV factors if both oil and gas active
530  const ADB& rv_perfcells = subset(state.rv, well_cells);
531  const ADB& rs_perfcells = subset(state.rs, well_cells);
532  const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
533 
534  const int oilpos = pu.phase_pos[Oil];
535  const int gaspos = pu.phase_pos[Gas];
536 
537  const ADB tmp_oil = (cmix_s[oilpos] - rv_perfcells * cmix_s[gaspos]) / d;
538  volumeRatio += tmp_oil / b_perfcells[oilpos];
539 
540  const ADB tmp_gas = (cmix_s[gaspos] - rs_perfcells * cmix_s[oilpos]) / d;
541  volumeRatio += tmp_gas / b_perfcells[gaspos];
542  }
543  else {
544  if ((*active_)[Oil]) {
545  const int oilpos = pu.phase_pos[Oil];
546  volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
547  }
548  if ((*active_)[Gas]) {
549  const int gaspos = pu.phase_pos[Gas];
550  volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
551  }
552  }
553 
554 
555  // injecting connections total volumerates at standard conditions
556  ADB cqt_is = cqt_i/volumeRatio;
557 
558  // connection phase volumerates at standard conditions
559  cq_s.resize(np, ADB::null());
560  for (int phase = 0; phase < np; ++phase) {
561  cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
562  }
563 
564  // check for dead wells (used in the well controll equations)
565  aliveWells = Vector::Constant(nw, 1.0);
566  for (int w = 0; w < nw; ++w) {
567  if (wbqt.value()[w] == 0) {
568  aliveWells[w] = 0.0;
569  }
570  }
571  }
572 
573 
574 
575 
576 
577  template <class SolutionState, class WellState>
578  void
579  StandardWells::
580  updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
581  const SolutionState& state,
582  WellState& xw) const
583  {
584  if ( !localWellsActive() )
585  {
586  // If there are no wells in the subdomain of the proces then
587  // cq_s has zero size and will cause a segmentation fault below.
588  return;
589  }
590 
591  // Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
592  const int np = wells().number_of_phases;
593  const int nw = wells().number_of_wells;
594  const int nperf = wells().well_connpos[nw];
595  Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
596  for (int phase = 1; phase < np; ++phase) {
597  cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
598  }
599  xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
600 
601  // Update the perforation pressures.
602  const Vector& cdp = wellPerforationPressureDiffs();
603  const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp;
604  xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
605  }
606 
607 
608 
609 
610 
611  template <class WellState>
612  void
613  StandardWells::
614  updateWellState(const Vector& dwells,
615  const double dpmaxrel,
616  WellState& well_state)
617  {
618  if( localWellsActive() )
619  {
620  const int np = wells().number_of_phases;
621  const int nw = wells().number_of_wells;
622 
623  // Extract parts of dwells corresponding to each part.
624  int varstart = 0;
625  const Vector dqs = subset(dwells, Span(np*nw, 1, varstart));
626  varstart += dqs.size();
627  const Vector dbhp = subset(dwells, Span(nw, 1, varstart));
628  varstart += dbhp.size();
629  assert(varstart == dwells.size());
630 
631 
632  // Qs update.
633  // Since we need to update the wellrates, that are ordered by wells,
634  // from dqs which are ordered by phase, the simplest is to compute
635  // dwr, which is the data from dqs but ordered by wells.
636  const DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
637  const Vector dwr = Eigen::Map<const Vector>(wwr.data(), nw*np);
638  const Vector wr_old = Eigen::Map<const Vector>(&well_state.wellRates()[0], nw*np);
639  const Vector wr = wr_old - dwr;
640  std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
641 
642  // Bhp update.
643  const Vector bhp_old = Eigen::Map<const Vector>(&well_state.bhp()[0], nw, 1);
644  const Vector dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel);
645  const Vector bhp = bhp_old - dbhp_limited;
646  std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
647 
648 
649  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
650  //Loop over all wells
651 #pragma omp parallel for schedule(static)
652  for (int w = 0; w < nw; ++w) {
653  const WellControls* wc = wells().ctrls[w];
654  const int nwc = well_controls_get_num(wc);
655  //Loop over all controls until we find a THP control
656  //that specifies what we need...
657  //Will only update THP for wells with THP control
658  for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
659  if (well_controls_iget_type(wc, ctrl_index) == THP) {
660  double aqua = 0.0;
661  double liquid = 0.0;
662  double vapour = 0.0;
663 
664  if ((*active_)[ Water ]) {
665  aqua = wr[w*np + pu.phase_pos[ Water ] ];
666  }
667  if ((*active_)[ Oil ]) {
668  liquid = wr[w*np + pu.phase_pos[ Oil ] ];
669  }
670  if ((*active_)[ Gas ]) {
671  vapour = wr[w*np + pu.phase_pos[ Gas ] ];
672  }
673 
674  double alq = well_controls_iget_alq(wc, ctrl_index);
675  int table_id = well_controls_iget_vfp(wc, ctrl_index);
676 
677  const WellType& well_type = wells().type[w];
678  const int perf = wells().well_connpos[w]; //first perforation.
679  if (well_type == INJECTOR) {
680  double dp = wellhelpers::computeHydrostaticCorrection(
681  wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
682  wellPerforationDensities()[perf], gravity_);
683 
684  well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
685  }
686  else if (well_type == PRODUCER) {
687  double dp = wellhelpers::computeHydrostaticCorrection(
688  wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
689  wellPerforationDensities()[perf], gravity_);
690 
691  well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
692  }
693  else {
694  OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
695  }
696 
697  //Assume only one THP control specified for each well
698  break;
699  }
700  }
701  }
702  }
703  }
704 
705 
706 
707 
708 
709  template <class WellState>
710  void
711  StandardWells::
712  updateWellControls(WellState& xw) const
713  {
714  wellhelpers::WellSwitchingLogger logger;
715 
716  if( !localWellsActive() ) return ;
717 
718  // Find, for each well, if any constraints are broken. If so,
719  // switch control to first broken constraint.
720  const int np = wells().number_of_phases;
721  const int nw = wells().number_of_wells;
722 #pragma omp parallel for schedule(dynamic)
723  for (int w = 0; w < nw; ++w) {
724  const WellControls* wc = wells().ctrls[w];
725  // The current control in the well state overrides
726  // the current control set in the Wells struct, which
727  // is instead treated as a default.
728  int current = xw.currentControls()[w];
729 
730  // Loop over all controls except the current one, and also
731  // skip any RESERVOIR_RATE controls, since we cannot
732  // handle those.
733  const int nwc = well_controls_get_num(wc);
734 
735  // There should be at least one control
736  assert(nwc != 0);
737 
738  bool constraint_violated = false;
739  int number_iterations = 0;
740  const int max_iterations = 2 * nwc; // maximum allowed iterations
741  do {
742  updateWellStateWithTarget(wc, current, w, xw);
743  int ctrl_index = 0;
744  for (; ctrl_index < nwc; ++ctrl_index) {
745  if (ctrl_index == current) {
746  // This is the currently used control, so it is
747  // used as an equation. So this is not used as an
748  // inequality constraint, and therefore skipped.
749  continue;
750  }
751  if (wellhelpers::constraintBroken(
752  xw.bhp(), xw.thp(), xw.wellRates(),
753  w, np, wells().type[w], wc, ctrl_index)) {
754  // ctrl_index will be the index of the broken constraint after the loop.
755  break;
756  }
757  }
758 
759 
760  if (ctrl_index != nwc) {
761  // Constraint number ctrl_index was broken, switch to it.
762  // We disregard terminal_ouput here as with it only messages
763  // for wells on one process will be printed.
764  logger.wellSwitched(wells().name[w],
765  well_controls_iget_type(wc, current),
766  well_controls_iget_type(wc, ctrl_index));
767 
768  xw.currentControls()[w] = ctrl_index;
769  current = xw.currentControls()[w];
770  constraint_violated = true;
771  } else {
772  constraint_violated = false;
773  }
774  ++number_iterations;
775 
776  if (number_iterations > max_iterations) {
777  OPM_THROW(Opm::NumericalProblem, "Could not find proper control within " << number_iterations << " iterations!");
778  break;
779  }
780  } while (constraint_violated);
781 
782 
783  if (wellCollection()->groupControlActive()) {
784 
785  // get well node in the well collection
786  WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
787 
788  // update whehter the well is under group control or individual control
789  if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
790  // under group control
791  well_node.setIndividualControl(false);
792  } else {
793  // individual control
794  well_node.setIndividualControl(true);
795  }
796  }
797  }
798 
799  }
800 
801 
802 
803 
804 
805  template <class SolutionState>
806  void
807  StandardWells::
808  addWellFluxEq(const std::vector<ADB>& cq_s,
809  const SolutionState& state,
810  LinearisedBlackoilResidual& residual)
811  {
812  if( !localWellsActive() )
813  {
814  // If there are no wells in the subdomain of the proces then
815  // cq_s has zero size and will cause a segmentation fault below.
816  return;
817  }
818 
819  const int np = wells().number_of_phases;
820  const int nw = wells().number_of_wells;
821  ADB qs = state.qs;
822  for (int phase = 0; phase < np; ++phase) {
823  qs -= superset(wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
824 
825  }
826 
827  residual.well_flux_eq = qs;
828  }
829 
830 
831 
832 
833 
834  template <class SolutionState, class WellState>
835  void
836  StandardWells::addWellControlEq(const SolutionState& state,
837  const WellState& xw,
838  const Vector& aliveWells,
839  LinearisedBlackoilResidual& residual)
840  {
841  if( ! localWellsActive() ) return;
842 
843  const int np = wells().number_of_phases;
844  const int nw = wells().number_of_wells;
845 
846  ADB aqua = ADB::constant(Vector::Zero(nw));
847  ADB liquid = ADB::constant(Vector::Zero(nw));
848  ADB vapour = ADB::constant(Vector::Zero(nw));
849 
850  if ((*active_)[Water]) {
851  aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw));
852  }
853  if ((*active_)[Oil]) {
854  liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw));
855  }
856  if ((*active_)[Gas]) {
857  vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw));
858  }
859 
860  //THP calculation variables
861  std::vector<int> inj_table_id(nw, -1);
862  std::vector<int> prod_table_id(nw, -1);
863  Vector thp_inj_target_v = Vector::Zero(nw);
864  Vector thp_prod_target_v = Vector::Zero(nw);
865  Vector alq_v = Vector::Zero(nw);
866 
867  //Hydrostatic correction variables
868  Vector rho_v = Vector::Zero(nw);
869  Vector vfp_ref_depth_v = Vector::Zero(nw);
870 
871  //Target vars
872  Vector bhp_targets = Vector::Zero(nw);
873  Vector rate_targets = Vector::Zero(nw);
874  Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
875 
876  //Selection variables
877  std::vector<int> bhp_elems;
878  std::vector<int> thp_inj_elems;
879  std::vector<int> thp_prod_elems;
880  std::vector<int> rate_elems;
881 
882  //Run through all wells to calculate BHP/RATE targets
883  //and gather info about current control
884  for (int w = 0; w < nw; ++w) {
885  auto wc = wells().ctrls[w];
886 
887  // The current control in the well state overrides
888  // the current control set in the Wells struct, which
889  // is instead treated as a default.
890  const int current = xw.currentControls()[w];
891 
892  switch (well_controls_iget_type(wc, current)) {
893  case BHP:
894  {
895  bhp_elems.push_back(w);
896  bhp_targets(w) = well_controls_iget_target(wc, current);
897  rate_targets(w) = -1e100;
898  }
899  break;
900 
901  case THP:
902  {
903  const int perf = wells().well_connpos[w];
904  rho_v[w] = wellPerforationDensities()[perf];
905 
906  const int table_id = well_controls_iget_vfp(wc, current);
907  const double target = well_controls_iget_target(wc, current);
908 
909  const WellType& well_type = wells().type[w];
910  if (well_type == INJECTOR) {
911  inj_table_id[w] = table_id;
912  thp_inj_target_v[w] = target;
913  alq_v[w] = -1e100;
914 
915  vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
916 
917  thp_inj_elems.push_back(w);
918  }
919  else if (well_type == PRODUCER) {
920  prod_table_id[w] = table_id;
921  thp_prod_target_v[w] = target;
922  alq_v[w] = well_controls_iget_alq(wc, current);
923 
924  vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
925 
926  thp_prod_elems.push_back(w);
927  }
928  else {
929  OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well");
930  }
931  bhp_targets(w) = -1e100;
932  rate_targets(w) = -1e100;
933  }
934  break;
935 
936  case RESERVOIR_RATE: // Intentional fall-through
937  case SURFACE_RATE:
938  {
939  rate_elems.push_back(w);
940  // RESERVOIR and SURFACE rates look the same, from a
941  // high-level point of view, in the system of
942  // simultaneous linear equations.
943 
944  const double* const distr =
945  well_controls_iget_distr(wc, current);
946 
947  for (int p = 0; p < np; ++p) {
948  rate_distr.insert(w, p*nw + w) = distr[p];
949  }
950 
951  bhp_targets(w) = -1.0e100;
952  rate_targets(w) = well_controls_iget_target(wc, current);
953  }
954  break;
955  }
956  }
957 
958  //Calculate BHP target from THP
959  const ADB thp_inj_target = ADB::constant(thp_inj_target_v);
960  const ADB thp_prod_target = ADB::constant(thp_prod_target_v);
961  const ADB alq = ADB::constant(alq_v);
962  const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target);
963  const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq);
964 
965  //Perform hydrostatic correction to computed targets
966  const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), gravity_);
967  const ADB dp = ADB::constant(dp_v);
968  const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw);
969  const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
970 
971  //Calculate residuals
972  const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj;
973  const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod;
974  const ADB bhp_residual = state.bhp - bhp_targets;
975  const ADB rate_residual = rate_distr * state.qs - rate_targets;
976 
977  //Select the right residual for each well
978  residual.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) +
979  superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) +
980  superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) +
981  superset(subset(rate_residual, rate_elems), rate_elems, nw);
982 
983  // For wells that are dead (not flowing), and therefore not communicating
984  // with the reservoir, we set the equation to be equal to the well's total
985  // flow. This will be a solution only if the target rate is also zero.
986  Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
987  for (int w = 0; w < nw; ++w) {
988  for (int phase = 0; phase < np; ++phase) {
989  rate_summer.insert(w, phase*nw + w) = 1.0;
990  }
991  }
992  Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
993  residual.well_eq = alive_selector.select(residual.well_eq, rate_summer * state.qs);
994  // OPM_AD_DUMP(residual_.well_eq);
995  }
996 
997 
998 
999 
1000 
1001  template <class SolutionState, class WellState>
1002  void
1003  StandardWells::computeWellPotentials(const std::vector<ADB>& mob_perfcells,
1004  const std::vector<ADB>& b_perfcells,
1005  const WellState& well_state,
1006  SolutionState& state0,
1007  std::vector<double>& well_potentials) const
1008  {
1009  const int nw = wells().number_of_wells;
1010  const int np = wells().number_of_phases;
1011  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
1012 
1013  Vector bhps = Vector::Zero(nw);
1014  for (int w = 0; w < nw; ++w) {
1015  const WellControls* ctrl = wells().ctrls[w];
1016  const int nwc = well_controls_get_num(ctrl);
1017  //Loop over all controls until we find a BHP control
1018  //or a THP control that specifies what we need.
1019  //Pick the value that gives the most restrictive flow
1020  for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
1021 
1022  if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
1023  bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
1024  }
1025 
1026  if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
1027  double aqua = 0.0;
1028  double liquid = 0.0;
1029  double vapour = 0.0;
1030 
1031  if ((*active_)[ Water ]) {
1032  aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
1033  }
1034  if ((*active_)[ Oil ]) {
1035  liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
1036  }
1037  if ((*active_)[ Gas ]) {
1038  vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
1039  }
1040 
1041  const int vfp = well_controls_iget_vfp(ctrl, ctrl_index);
1042  const double& thp = well_controls_iget_target(ctrl, ctrl_index);
1043  const double& alq = well_controls_iget_alq(ctrl, ctrl_index);
1044 
1045  //Set *BHP* target by calculating bhp from THP
1046  const WellType& well_type = wells().type[w];
1047  const int perf = wells().well_connpos[w]; //first perforation
1048 
1049  if (well_type == INJECTOR) {
1050  double dp = wellhelpers::computeHydrostaticCorrection(
1051  wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
1052  wellPerforationDensities()[perf], gravity_);
1053  const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
1054  // apply the strictest of the bhp controlls i.e. smallest bhp for injectors
1055  if ( bhp < bhps[w]) {
1056  bhps[w] = bhp;
1057  }
1058  }
1059  else if (well_type == PRODUCER) {
1060  double dp = wellhelpers::computeHydrostaticCorrection(
1061  wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
1062  wellPerforationDensities()[perf], gravity_);
1063 
1064  const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
1065  // apply the strictest of the bhp controlls i.e. largest bhp for producers
1066  if ( bhp > bhps[w]) {
1067  bhps[w] = bhp;
1068  }
1069  }
1070  else {
1071  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
1072  }
1073  }
1074  }
1075 
1076  }
1077 
1078  // use bhp limit from control
1079  state0.bhp = ADB::constant(bhps);
1080 
1081  // compute well potentials
1082  Vector aliveWells;
1083  std::vector<ADB> perf_potentials;
1084  computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, perf_potentials);
1085 
1086  well_potentials.resize(nw * np, 0.0);
1087 
1088  for (int p = 0; p < np; ++p) {
1089  for (int w = 0; w < nw; ++w) {
1090  for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) {
1091  well_potentials[w * np + p] += perf_potentials[p].value()[perf];
1092  }
1093  }
1094  }
1095  }
1096 
1097 
1098 
1099 
1100 
1101  void
1102  StandardWells::variableStateWellIndices(std::vector<int>& indices,
1103  int& next) const
1104  {
1105  indices[Qs] = next++;
1106  indices[Bhp] = next++;
1107  }
1108 
1109 
1110 
1111 
1112 
1113  template <class SolutionState>
1114  void
1115  StandardWells::
1116  variableStateExtractWellsVars(const std::vector<int>& indices,
1117  std::vector<ADB>& vars,
1118  SolutionState& state) const
1119  {
1120  // Qs.
1121  state.qs = std::move(vars[indices[Qs]]);
1122 
1123  // Bhp.
1124  state.bhp = std::move(vars[indices[Bhp]]);
1125  }
1126 
1127 
1128 
1129 
1130 
1131  std::vector<int>
1132  StandardWells::variableWellStateIndices() const
1133  {
1134  // Black oil model standard is 5 equation.
1135  // For the pure well solve, only the well equations are picked.
1136  std::vector<int> indices(5, -1);
1137  int next = 0;
1138 
1139  variableStateWellIndices(indices, next);
1140 
1141  assert(next == 2);
1142  return indices;
1143  }
1144 
1145 
1146 
1147 
1148 
1149  template <class WellState>
1150  void
1151  StandardWells::variableWellStateInitials(const WellState& xw,
1152  std::vector<Vector>& vars0) const
1153  {
1154  // Initial well rates.
1155  if ( localWellsActive() )
1156  {
1157  // Need to reshuffle well rates, from phase running fastest
1158  // to wells running fastest.
1159  const int nw = wells().number_of_wells;
1160  const int np = wells().number_of_phases;
1161 
1162  // The transpose() below switches the ordering.
1163  const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
1164  const Vector qs = Eigen::Map<const V>(wrates.data(), nw*np);
1165  vars0.push_back(qs);
1166 
1167  // Initial well bottom-hole pressure.
1168  assert (not xw.bhp().empty());
1169  const Vector bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
1170  vars0.push_back(bhp);
1171  }
1172  else
1173  {
1174  // push null states for qs and bhp
1175  vars0.push_back(V());
1176  vars0.push_back(V());
1177  }
1178  }
1179 
1180 
1181 
1182 
1183 
1184  void
1186  {
1187  store_well_perforation_fluxes_ = store_fluxes;
1188  }
1189 
1190 
1191 
1192 
1193 
1194  const StandardWells::Vector&
1196  {
1197  assert(store_well_perforation_fluxes_);
1198  return well_perforation_fluxes_;
1199  }
1200 
1201 
1202 
1203 
1204 
1205 
1206  template<class WellState>
1207  void
1209  updateListEconLimited(const Schedule& schedule,
1210  const int current_step,
1211  const Wells* wells_struct,
1212  const WellState& well_state,
1213  DynamicListEconLimited& list_econ_limited) const
1214  {
1215  // wells_struct may be null pointer if there are no wells in process domain
1216  const int nw = ( wells_struct ) ? wells_struct->number_of_wells : 0;
1217 
1218  for (int w = 0; w < nw; ++w) {
1219  // flag to check if the mim oil/gas rate limit is violated
1220  bool rate_limit_violated = false;
1221  const std::string& well_name = wells_struct->name[w];
1222  const Well* well_ecl = schedule.getWell(well_name);
1223  const WellEconProductionLimits& econ_production_limits = well_ecl->getEconProductionLimits(current_step);
1224 
1225  // economic limits only apply for production wells.
1226  if (wells_struct->type[w] != PRODUCER) {
1227  continue;
1228  }
1229 
1230  // if no limit is effective here, then continue to the next well
1231  if ( !econ_production_limits.onAnyEffectiveLimit() ) {
1232  continue;
1233  }
1234  // for the moment, we only handle rate limits, not handling potential limits
1235  // the potential limits should not be difficult to add
1236  const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit();
1237  if (quantity_limit == WellEcon::POTN) {
1238  const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n")
1239  + std::string("All the limits will be evaluated based on RATE. ");
1240  OpmLog::warning("NOT_SUPPORTING_POTN", msg);
1241  }
1242 
1243  const WellMapType& well_map = well_state.wellMap();
1244  const typename WellMapType::const_iterator i_well = well_map.find(well_name);
1245  assert(i_well != well_map.end()); // should always be found?
1246  const WellMapEntryType& map_entry = i_well->second;
1247  const int well_number = map_entry[0];
1248 
1249  if (econ_production_limits.onAnyRateLimit()) {
1250  rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state, well_number);
1251  }
1252 
1253  if (rate_limit_violated) {
1254  if (econ_production_limits.endRun()) {
1255  const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n")
1256  + std::string("the program will keep running after ") + well_name + std::string(" is closed");
1257  OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message);
1258  }
1259 
1260  if (econ_production_limits.validFollowonWell()) {
1261  OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
1262  }
1263 
1264  if (well_ecl->getAutomaticShutIn()) {
1265  list_econ_limited.addShutWell(well_name);
1266  const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to economic limit");
1267  OpmLog::info(msg);
1268  } else {
1269  list_econ_limited.addStoppedWell(well_name);
1270  const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to economic limit");
1271  OpmLog::info(msg);
1272  }
1273  // the well is closed, not need to check other limits
1274  continue;
1275  }
1276 
1277  // checking for ratio related limits, mostly all kinds of ratio.
1278  bool ratio_limits_violated = false;
1279  RatioCheckTuple ratio_check_return;
1280 
1281  if (econ_production_limits.onAnyRatioLimit()) {
1282  ratio_check_return = checkRatioEconLimits(econ_production_limits, well_state, map_entry);
1283  ratio_limits_violated = std::get<0>(ratio_check_return);
1284  }
1285 
1286  if (ratio_limits_violated) {
1287  const bool last_connection = std::get<1>(ratio_check_return);
1288  const int worst_offending_connection = std::get<2>(ratio_check_return);
1289 
1290  const int perf_start = map_entry[1];
1291 
1292  assert((worst_offending_connection >= 0) && (worst_offending_connection < map_entry[2]));
1293 
1294  const int cell_worst_offending_connection = wells_struct->well_cells[perf_start + worst_offending_connection];
1295  list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection);
1296  const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ")
1297  + well_name + std::string(" will be closed due to economic limit");
1298  OpmLog::info(msg);
1299 
1300  if (last_connection) {
1301  list_econ_limited.addShutWell(well_name);
1302  const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed");
1303  OpmLog::info(msg2);
1304  }
1305  }
1306 
1307  }
1308  }
1309 
1310 
1311 
1312 
1313 
1314  template <class WellState>
1315  bool
1316  StandardWells::
1317  checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
1318  const WellState& well_state,
1319  const int well_number) const
1320  {
1321  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
1322  const int np = well_state.numPhases();
1323 
1324  if (econ_production_limits.onMinOilRate()) {
1325  assert((*active_)[Oil]);
1326  const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
1327  const double min_oil_rate = econ_production_limits.minOilRate();
1328  if (std::abs(oil_rate) < min_oil_rate) {
1329  return true;
1330  }
1331  }
1332 
1333  if (econ_production_limits.onMinGasRate() ) {
1334  assert((*active_)[Gas]);
1335  const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ];
1336  const double min_gas_rate = econ_production_limits.minGasRate();
1337  if (std::abs(gas_rate) < min_gas_rate) {
1338  return true;
1339  }
1340  }
1341 
1342  if (econ_production_limits.onMinLiquidRate() ) {
1343  assert((*active_)[Oil]);
1344  assert((*active_)[Water]);
1345  const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
1346  const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
1347  const double liquid_rate = oil_rate + water_rate;
1348  const double min_liquid_rate = econ_production_limits.minLiquidRate();
1349  if (std::abs(liquid_rate) < min_liquid_rate) {
1350  return true;
1351  }
1352  }
1353 
1354  if (econ_production_limits.onMinReservoirFluidRate()) {
1355  OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
1356  }
1357 
1358  return false;
1359  }
1360 
1361 
1362 
1363 
1364 
1365  template <class WellState>
1366  StandardWells::RatioCheckTuple
1367  StandardWells::
1368  checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
1369  const WellState& well_state,
1370  const WellMapEntryType& map_entry) const
1371  {
1372  // TODO: not sure how to define the worst-offending connection when more than one
1373  // ratio related limit is violated.
1374  // The defintion used here is that we define the violation extent based on the
1375  // ratio between the value and the corresponding limit.
1376  // For each violated limit, we decide the worst-offending connection separately.
1377  // Among the worst-offending connections, we use the one has the biggest violation
1378  // extent.
1379 
1380  bool any_limit_violated = false;
1381  bool last_connection = false;
1382  int worst_offending_connection = INVALIDCONNECTION;
1383  double violation_extent = -1.0;
1384 
1385  if (econ_production_limits.onMaxWaterCut()) {
1386  const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry);
1387  bool water_cut_violated = std::get<0>(water_cut_return);
1388  if (water_cut_violated) {
1389  any_limit_violated = true;
1390  const double violation_extent_water_cut = std::get<3>(water_cut_return);
1391  if (violation_extent_water_cut > violation_extent) {
1392  violation_extent = violation_extent_water_cut;
1393  worst_offending_connection = std::get<2>(water_cut_return);
1394  last_connection = std::get<1>(water_cut_return);
1395  }
1396  }
1397  }
1398 
1399  if (econ_production_limits.onMaxGasOilRatio()) {
1400  OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
1401  }
1402 
1403  if (econ_production_limits.onMaxWaterGasRatio()) {
1404  OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
1405  }
1406 
1407  if (econ_production_limits.onMaxGasLiquidRatio()) {
1408  OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
1409  }
1410 
1411  if (any_limit_violated) {
1412  assert(worst_offending_connection >=0);
1413  assert(violation_extent > 1.);
1414  }
1415 
1416  return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
1417  }
1418 
1419 
1420 
1421 
1422 
1423  template <class WellState>
1424  StandardWells::RatioCheckTuple
1425  StandardWells::
1426  checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
1427  const WellState& well_state,
1428  const WellMapEntryType& map_entry) const
1429  {
1430  bool water_cut_limit_violated = false;
1431  int worst_offending_connection = INVALIDCONNECTION;
1432  bool last_connection = false;
1433  double violation_extent = -1.0;
1434 
1435  const int np = well_state.numPhases();
1436  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
1437  const int well_number = map_entry[0];
1438 
1439  assert((*active_)[Oil]);
1440  assert((*active_)[Water]);
1441 
1442  const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
1443  const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
1444  const double liquid_rate = oil_rate + water_rate;
1445  double water_cut;
1446  if (std::abs(liquid_rate) != 0.) {
1447  water_cut = water_rate / liquid_rate;
1448  } else {
1449  water_cut = 0.0;
1450  }
1451 
1452  const double max_water_cut_limit = econ_production_limits.maxWaterCut();
1453  if (water_cut > max_water_cut_limit) {
1454  water_cut_limit_violated = true;
1455  }
1456 
1457  if (water_cut_limit_violated) {
1458  // need to handle the worst_offending_connection
1459  const int perf_start = map_entry[1];
1460  const int perf_number = map_entry[2];
1461 
1462  std::vector<double> water_cut_perf(perf_number);
1463  for (int perf = 0; perf < perf_number; ++perf) {
1464  const int i_perf = perf_start + perf;
1465  const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
1466  const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
1467  const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
1468  if (std::abs(liquid_perf_rate) != 0.) {
1469  water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
1470  } else {
1471  water_cut_perf[perf] = 0.;
1472  }
1473  }
1474 
1475  last_connection = (perf_number == 1);
1476  if (last_connection) {
1477  worst_offending_connection = 0;
1478  violation_extent = water_cut_perf[0] / max_water_cut_limit;
1479  return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
1480  }
1481 
1482  double max_water_cut_perf = 0.;
1483  for (int perf = 0; perf < perf_number; ++perf) {
1484  if (water_cut_perf[perf] > max_water_cut_perf) {
1485  worst_offending_connection = perf;
1486  max_water_cut_perf = water_cut_perf[perf];
1487  }
1488  }
1489 
1490  assert(max_water_cut_perf != 0.);
1491  assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
1492 
1493  violation_extent = max_water_cut_perf / max_water_cut_limit;
1494  }
1495 
1496  return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
1497  }
1498 
1499 
1500 
1501 
1502 
1503  WellCollection* StandardWells::wellCollection() const
1504  {
1505  return well_collection_;
1506  }
1507 
1508 
1509 
1510 
1511 
1512  void StandardWells::calculateEfficiencyFactors()
1513  {
1514  if ( !localWellsActive() ) {
1515  return;
1516  }
1517  // get efficiency factor for each well first
1518  const int nw = wells_->number_of_wells;
1519 
1520  Vector well_efficiency_factors = Vector::Ones(nw);
1521 
1522  for (int w = 0; w < nw; ++w) {
1523  const std::string well_name = wells_->name[w];
1524  const WellNode& well_node = well_collection_->findWellNode(well_name);
1525 
1526  well_efficiency_factors(w) = well_node.getAccumulativeEfficiencyFactor();
1527  }
1528 
1529  // map them to the perforation.
1530  well_perforation_efficiency_factors_ = wellOps().w2p * well_efficiency_factors.matrix();
1531  }
1532 
1533 
1534 
1535 
1536 
1537  const StandardWells::Vector&
1538  StandardWells::wellPerfEfficiencyFactors() const
1539  {
1540  return well_perforation_efficiency_factors_;
1541  }
1542 
1543 
1544 
1545 
1546 
1547  template <class WellState>
1548  void
1549  StandardWells::
1550  updateWellStateWithTarget(const WellControls* wc,
1551  const int current,
1552  const int well_index,
1553  WellState& xw) const
1554  {
1555  const int np = wells().number_of_phases;
1556 
1557  // Updating well state and primary variables.
1558  // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
1559  const double target = well_controls_iget_target(wc, current);
1560  const double* distr = well_controls_iget_distr(wc, current);
1561  switch (well_controls_iget_type(wc, current)) {
1562  case BHP:
1563  xw.bhp()[well_index] = target;
1564  break;
1565  case THP: {
1566  double aqua = 0.0;
1567  double liquid = 0.0;
1568  double vapour = 0.0;
1569 
1570  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
1571 
1572  if ((*active_)[ Water ]) {
1573  aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
1574  }
1575  if ((*active_)[ Oil ]) {
1576  liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
1577  }
1578  if ((*active_)[ Gas ]) {
1579  vapour = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
1580  }
1581 
1582  const int vfp = well_controls_iget_vfp(wc, current);
1583  const double& thp = well_controls_iget_target(wc, current);
1584  const double& alq = well_controls_iget_alq(wc, current);
1585 
1586  //Set *BHP* target by calculating bhp from THP
1587  const WellType& well_type = wells().type[well_index];
1588  // pick the density in the top layer
1589  const int perf = wells().well_connpos[well_index];
1590  const double rho = well_perforation_densities_[perf];
1591 
1592  if (well_type == INJECTOR) {
1593  double dp = wellhelpers::computeHydrostaticCorrection(
1594  wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
1595  rho, gravity_);
1596 
1597  xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
1598  }
1599  else if (well_type == PRODUCER) {
1600  double dp = wellhelpers::computeHydrostaticCorrection(
1601  wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
1602  rho, gravity_);
1603 
1604  xw.bhp()[well_index] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
1605  }
1606  else {
1607  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
1608  }
1609  break;
1610  }
1611 
1612  case RESERVOIR_RATE:
1613  // No direct change to any observable quantity at
1614  // surface condition. In this case, use existing
1615  // flow rates as initial conditions as reservoir
1616  // rate acts only in aggregate.
1617  break;
1618 
1619  case SURFACE_RATE:
1620  // assign target value as initial guess for injectors and
1621  // single phase producers (orat, grat, wrat)
1622  const WellType& well_type = wells().type[well_index];
1623  if (well_type == INJECTOR) {
1624  for (int phase = 0; phase < np; ++phase) {
1625  const double& compi = wells().comp_frac[np * well_index + phase];
1626  if (compi > 0.0) {
1627  xw.wellRates()[np * well_index + phase] = target * compi;
1628  }
1629  }
1630  } else if (well_type == PRODUCER) {
1631 
1632  // only set target as initial rates for single phase
1633  // producers. (orat, grat and wrat, and not lrat)
1634  // lrat will result in numPhasesWithTargetsUnderThisControl == 2
1635  int numPhasesWithTargetsUnderThisControl = 0;
1636  for (int phase = 0; phase < np; ++phase) {
1637  if (distr[phase] > 0.0) {
1638  numPhasesWithTargetsUnderThisControl += 1;
1639  }
1640  }
1641  for (int phase = 0; phase < np; ++phase) {
1642  if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
1643  xw.wellRates()[np * well_index + phase] = target * distr[phase];
1644  }
1645  }
1646  } else {
1647  OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
1648  }
1649 
1650  break;
1651  }
1652  }
1653 
1654 
1655 } // namespace Opm
ADB bWat(const ADB &pw, const ADB &T, const Cells &cells) const
Water formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:446
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
V rsSat(const V &po, const Cells &cells) const
Bubble point curve for Rs as function of oil pressure.
Vector & wellPerforationDensities()
Density of each well perforation.
Definition: StandardWells_impl.hpp:184
ADB bhp(const std::vector< int > &table_id, const Wells &wells, const ADB &qs, const ADB &thp) const
Linear interpolation of bhp as function of the input parameters.
Definition: VFPInjProperties.cpp:71
ADB rvSat(const ADB &po, const Cells &cells) const
Condensation curve for Rv as function of oil pressure.
Definition: BlackoilPropsAdFromDeck.cpp:688
double thp(int table_id, const double &aqua, const double &liquid, const double &vapour, const double &bhp, const double &alq) const
Linear interpolation of thp as a function of the input parameters.
Definition: VFPProdProperties.cpp:192
bool localWellsActive() const
return true if wells are available on this process
Definition: StandardWells_impl.hpp:148
bool wellsActive() const
return true if wells are available in the reservoir
Definition: StandardWells_impl.hpp:130
const VFPInjTable * getTable(const int table_id) const
Returns the table associated with the ID, or throws an exception if the table does not exist...
Definition: VFPInjProperties.cpp:218
Vector & wellPerforationPressureDiffs()
Diff to bhp for each well perforation.
Definition: StandardWells_impl.hpp:204
void setStoreWellPerforationFluxesFlag(const bool store_fluxes)
If set, computeWellFlux() will additionally store the total reservoir volume perforation fluxes...
Definition: StandardWells_impl.hpp:1185
double thp(int table_id, const double &aqua, const double &liquid, const double &vapour, const double &bhp) const
Linear interpolation of thp as a function of the input parameters.
Definition: VFPInjProperties.cpp:183
ADB bOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const
Oil formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:493
const VFPInjProperties * getInj() const
Returns the VFP properties for injection wells.
Definition: VFPProperties.hpp:61
const VFPProdTable * getTable(const int table_id) const
Returns the table associated with the ID, or throws an exception if the table does not exist...
Definition: VFPProdProperties.cpp:234
ADB bGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const
Gas formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:566
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
V surfaceDensity(const int phaseIdx, const Cells &cells) const
Densities of stock components at surface conditions.
Definition: BlackoilPropsAdFromDeck.cpp:242
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
static std::vector< double > computeConnectionPressureDelta(const Wells &wells, const std::vector< double > &z_perf, const std::vector< double > &dens_perf, const double gravity)
Compute pressure deltas.
Definition: WellDensitySegmented.cpp:271
static std::vector< double > computeConnectionDensities(const Wells &wells, const PhaseUsage &phase_usage, const std::vector< double > &perfComponentRates, const std::vector< double > &b_perf, const std::vector< double > &rsmax_perf, const std::vector< double > &rvmax_perf, const std::vector< double > &surf_dens_perf)
Compute well segment densities Notation: N = number of perforations, C = number of components...
Definition: WellDensitySegmented.cpp:32
ADB bhp(const std::vector< int > &table_id, const Wells &wells, const ADB &qs, const ADB &thp, const ADB &alq) const
Linear interpolation of bhp as function of the input parameters.
Definition: VFPProdProperties.cpp:58
const VFPProdProperties * getProd() const
Returns the VFP properties for production wells.
Definition: VFPProperties.hpp:68
PhaseUsage phaseUsage() const
Definition: BlackoilPropsAdFromDeck.cpp:231
void updateListEconLimited(const Schedule &schedule, const int current_step, const Wells *wells, const WellState &well_state, DynamicListEconLimited &list_econ_limited) const
upate the dynamic lists related to economic limits
Definition: StandardWells_impl.hpp:1209
const Vector & getStoredWellPerforationFluxes() const
Retrieves the stored fluxes.
Definition: StandardWells_impl.hpp:1195
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719