MultisegmentWells_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 #ifndef OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
22 #define OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
23 
24 
25 namespace Opm
26 {
27 
28 
29 
30  namespace wellhelpers {
31 
32  using ADB = MultisegmentWells::ADB;
33  using Vector = MultisegmentWells::Vector;
34 
35  inline
36  ADB onlyWellDerivs(const ADB& x)
37  {
38  Vector val = x.value();
39  const int nb = x.numBlocks();
40  if (nb < 2) {
41  OPM_THROW(std::logic_error, "Called onlyWellDerivs() with argument that has " << nb << " blocks.");
42  }
43  std::vector<ADB::M> derivs = { x.derivative()[nb - 2], x.derivative()[nb - 1] };
44  return ADB::function(std::move(val), std::move(derivs));
45  }
46  }
47 
48 
49 
50 
51 
52  template <class ReservoirResidualQuant, class SolutionState>
53  void
54  MultisegmentWells::
55  extractWellPerfProperties(const SolutionState& /* state */,
56  const std::vector<ReservoirResidualQuant>& rq,
57  std::vector<ADB>& mob_perfcells,
58  std::vector<ADB>& b_perfcells) const
59  {
60  // If we have wells, extract the mobilities and b-factors for
61  // the well-perforated cells.
62  if ( !localWellsActive() ) {
63  mob_perfcells.clear();
64  b_perfcells.clear();
65  return;
66  } else {
67  const std::vector<int>& well_cells = wellOps().well_cells;
68  mob_perfcells.resize(num_phases_, ADB::null());
69  b_perfcells.resize(num_phases_, ADB::null());
70  for (int phase = 0; phase < num_phases_; ++phase) {
71  mob_perfcells[phase] = subset(rq[phase].mob, well_cells);
72  b_perfcells[phase] = subset(rq[phase].b, well_cells);
73  }
74  }
75  }
76 
77 
78 
79 
80 
81  template <class WellState>
82  void
83  MultisegmentWells::
84  updateWellState(const Vector& dwells,
85  const double dpmaxrel,
86  WellState& well_state) const
87  {
88  if (!msWells().empty())
89  {
90  const int nw = msWells().size();
91  const int nseg_total = nseg_total_;
92  const int np = numPhases();
93 
94  // Extract parts of dwells corresponding to each part.
95  int varstart = 0;
96  const Vector dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart));
97  varstart += dsegqs.size();
98  const Vector dsegp = subset(dwells, Span(nseg_total, 1, varstart));
99  varstart += dsegp.size();
100  assert(varstart == dwells.size());
101 
102 
103  // segment phase rates update
104  // in dwells, the phase rates are ordered by phase.
105  // while in WellStateMultiSegment, the phase rates are ordered by segments
106  const DataBlock wsr = Eigen::Map<const DataBlock>(dsegqs.data(), np, nseg_total).transpose();
107  const Vector dwsr = Eigen::Map<const Vector>(wsr.data(), nseg_total * np);
108  const Vector wsr_old = Eigen::Map<const Vector>(&well_state.segPhaseRates()[0], nseg_total * np);
109  const Vector sr = wsr_old - dwsr;
110  std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin());
111 
112 
113  // segment pressure updates
114  const Vector segp_old = Eigen::Map<const Vector>(&well_state.segPress()[0], nseg_total, 1);
115  // TODO: applying the pressure change limiter to all the segments, not sure if it is the correct thing to do
116  const Vector dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel);
117  const Vector segp = segp_old - dsegp_limited;
118  std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin());
119 
120  // update the well rates and bhps, which are not anymore primary vabriables.
121  // they are updated directly from the updated segment phase rates and segment pressures.
122 
123  // Bhp update.
124  Vector bhp = Vector::Zero(nw);
125  Vector wr = Vector::Zero(nw * np);
126  // it is better to use subset
127 
128  int start_segment = 0;
129  for (int w = 0; w < nw; ++w) {
130  bhp[w] = well_state.segPress()[start_segment];
131  // insert can be faster
132  for (int p = 0; p < np; ++p) {
133  wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment];
134  }
135 
136  const int nseg = msWells()[w]->numberOfSegments();
137  start_segment += nseg;
138  }
139 
140  assert(start_segment == nseg_total);
141  std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
142  std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
143 
144  // TODO: handling the THP control related.
145  }
146  }
147 
148 
149 
150 
151 
152  template <class SolutionState>
153  void
154  MultisegmentWells::
155  computeWellFlux(const SolutionState& state,
156  const std::vector<ADB>& mob_perfcells,
157  const std::vector<ADB>& b_perfcells,
158  Vector& aliveWells,
159  std::vector<ADB>& cq_s) const
160  {
161  if (msWells().size() == 0) return;
162 
163  const int np = numPhases();
164  const int nw = msWells().size();
165 
166  aliveWells = Vector::Constant(nw, 1.0);
167 
168  const int nseg = nseg_total_;
169  const int nperf = nperf_total_;
170 
171  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
172 
173  cq_s.resize(np, ADB::null());
174 
175  {
176  const Vector& Tw = wellOps().conn_trans_factors;
177  const std::vector<int>& well_cells = wellOps().well_cells;
178 
179  // determining in-flow (towards well-bore) or out-flow (towards reservoir)
180  // for mutli-segmented wells and non-segmented wells, the calculation of the drawdown are different.
181  const ADB& p_perfcells = subset(state.pressure, well_cells);
182  const ADB& rs_perfcells = subset(state.rs, well_cells);
183  const ADB& rv_perfcells = subset(state.rv, well_cells);
184 
185  const ADB& seg_pressures = state.segp;
186 
187  const ADB seg_pressures_perf = wellOps().s2p * seg_pressures;
188 
189  // Create selector for perforations of multi-segment vs. regular wells.
190  Vector is_multisegment_well(nw);
191  for (int w = 0; w < nw; ++w) {
192  is_multisegment_well[w] = double(msWells()[w]->isMultiSegmented());
193  }
194  // Take one flag per well and expand to one flag per perforation.
195  Vector is_multisegment_perf = wellOps().w2p * is_multisegment_well.matrix();
196  Selector<double> msperf_selector(is_multisegment_perf, Selector<double>::NotEqualZero);
197 
198  // Compute drawdown.
199  ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_,
200  ADB::constant(well_perforation_pressure_diffs_));
201  const Vector h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, Vector::Zero(nperf));
202 
203  // Special handling for when we are called from solveWellEq().
204  // TODO: restructure to eliminate need for special treatmemt.
205  if ((h_nc.numBlocks() != 0) && (h_nc.numBlocks() != seg_pressures_perf.numBlocks())) {
206  assert(seg_pressures_perf.numBlocks() == 2);
207  assert(h_nc.numBlocks() > 2);
208  h_nc = wellhelpers::onlyWellDerivs(h_nc);
209  assert(h_nc.numBlocks() == 2);
210  }
211 
212  ADB drawdown = (p_perfcells + h_cj - seg_pressures_perf - h_nc);
213 
214  // selects injection perforations
215  Vector selectInjectingPerforations = Vector::Zero(nperf);
216  // selects producing perforations
217  Vector selectProducingPerforations = Vector::Zero(nperf);
218  for (int c = 0; c < nperf; ++c){
219  if (drawdown.value()[c] < 0)
220  selectInjectingPerforations[c] = 1;
221  else
222  selectProducingPerforations[c] = 1;
223  }
224 
225  // handling flow into wellbore
226  // maybe there are something to do there make the procedure easier.
227  std::vector<ADB> cq_ps(np, ADB::null());
228  for (int phase = 0; phase < np; ++phase) {
229  const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
230  cq_ps[phase] = b_perfcells[phase] * cq_p;
231  }
232 
233  if ((*active_)[Oil] && (*active_)[Gas]) {
234  const int oilpos = pu.phase_pos[Oil];
235  const int gaspos = pu.phase_pos[Gas];
236  const ADB cq_psOil = cq_ps[oilpos];
237  const ADB cq_psGas = cq_ps[gaspos];
238  cq_ps[gaspos] += rs_perfcells * cq_psOil;
239  cq_ps[oilpos] += rv_perfcells * cq_psGas;
240  }
241 
242  // hadling flow out from wellbore
243  ADB total_mob = mob_perfcells[0];
244  for (int phase = 1; phase < np; ++phase) {
245  total_mob += mob_perfcells[phase];
246  }
247 
248  // injection perforations total volume rates
249  const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
250 
251  // compute wellbore mixture for injecting perforations
252  // The wellbore mixture depends on the inflow from the reservoir
253  // and the well injection rates.
254  // TODO: should this based on the segments?
255  // TODO: for the usual wells, the well rates are the sum of the perforations.
256  // TODO: for multi-segmented wells, the segment rates are not the sum of the perforations.
257 
258  // TODO: two options here
259  // TODO: 1. for each segment, only the inflow from the perforations related to this segment are considered.
260  // TODO: 2. for each segment, the inflow from the perforrations related to this segment and also all the inflow
261  // TODO: from the upstreaming sgments and their perforations need to be considered.
262  // TODO: This way can be the more consistent way, while let us begin with the first option. The second option
263  // TODO: involves one operations that are not valid now. (i.e. how to transverse from the leaves to the root,
264  // TODO: although we can begin from the brutal force way)
265 
266  // TODO: stop using msWells() here.
267  std::vector<ADB> wbq(np, ADB::null());
268  ADB wbqt = ADB::constant(Vector::Zero(nseg));
269 
270  const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
271 
272  for (int phase = 0; phase < np; ++phase) {
273  const ADB& q_ps = wellOps().p2s * cq_ps[phase];
274  const ADB& q_s = subset(state.segqs, Span(nseg, 1, phase * nseg));
275  Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
276 
277  const int pos = pu.phase_pos[phase];
278 
279  // this is per segment
280  wbq[phase] = (wellOps().w2s * ADB::constant(compi.col(pos)) * injectingPhase_selector.select(q_s, ADB::constant(Vector::Zero(nseg)))) - q_ps;
281 
282  // TODO: it should be a single value for this certain well.
283  // TODO: it need to be changed later to handle things more consistently
284  // or there should be an earsier way to decide if the well is dead.
285  wbqt += wbq[phase];
286  }
287 
288  // Set aliveWells.
289  // the first value of the wbqt is the one to decide if the well is dead
290  // or there should be some dead segments?
291  {
292  int topseg = 0;
293  for (int w = 0; w < nw; ++w) {
294  if (wbqt.value()[topseg] == 0.0) { // yes we really mean == here, no fuzzyness
295  aliveWells[w] = 0.0;
296  }
297  topseg += msWells()[w]->numberOfSegments();
298  }
299  }
300 
301  // compute wellbore mixture at standard conditions.
302  // before, the determination of alive wells is based on wells.
303  // now, will there be any dead segment? I think no.
304  // TODO: it is not clear if the cmix_s should be based on segment or the well
305  std::vector<ADB> cmix_s(np, ADB::null());
306  Selector<double> aliveWells_selector(aliveWells, Selector<double>::NotEqualZero);
307  for (int phase = 0; phase < np; ++phase) {
308  const int pos = pu.phase_pos[phase];
309  const ADB phase_fraction = wellOps().topseg2w * (wbq[phase] / wbqt);
310  cmix_s[phase] = wellOps().w2p * aliveWells_selector.select(phase_fraction, ADB::constant(compi.col(pos)));
311  }
312 
313  // compute volume ration between connection at standard conditions
314  ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
315  const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
316 
317  for (int phase = 0; phase < np; ++phase) {
318  ADB tmp = cmix_s[phase];
319  if (phase == Oil && (*active_)[Gas]) {
320  const int gaspos = pu.phase_pos[Gas];
321  tmp = (tmp - rv_perfcells * cmix_s[gaspos]) / d;
322  }
323  if (phase == Gas && (*active_)[Oil]) {
324  const int oilpos = pu.phase_pos[Oil];
325  tmp = (tmp - rs_perfcells * cmix_s[oilpos]) / d;
326  }
327  volumeRatio += tmp / b_perfcells[phase];
328  }
329 
330  // injecting connections total volumerates at standard conditions
331  ADB cqt_is = cqt_i/volumeRatio;
332 
333  // connection phase volumerates at standard conditions
334  for (int phase = 0; phase < np; ++phase) {
335  cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
336  }
337  }
338  }
339 
340 
341 
342 
343 
344  template <class SolutionState, class WellState>
345  void
346  MultisegmentWells::
347  updatePerfPhaseRatesAndPressures(const std::vector<ADB>& cq_s,
348  const SolutionState& state,
349  WellState& xw) const
350  {
351  if ( !localWellsActive() ) {
352  return;
353  }
354  // Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
355  const int np = numPhases();
356  const int nw = numWells();
357 
358  Vector cq = superset(cq_s[0].value(), Span(nperf_total_, np, 0), nperf_total_ * np);
359  for (int phase = 1; phase < np; ++phase) {
360  cq += superset(cq_s[phase].value(), Span(nperf_total_, np, phase), nperf_total_ * np);
361  }
362  xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf_total_ * np);
363 
364  // Update the perforation pressures for usual wells first to recover the resutls
365  // without mutlti segment wells. For segment wells, it has not been decided if
366  // we need th concept of preforation pressures
367  xw.perfPress().resize(nperf_total_, -1.e100);
368 
369  const Vector& cdp = well_perforation_pressure_diffs_;
370  int start_segment = 0;
371  int start_perforation = 0;
372  for (int i = 0; i < nw; ++i) {
373  WellMultiSegmentConstPtr well = wells_multisegment_[i];
374  const int nperf = well->numberOfPerforations();
375  const int nseg = well->numberOfSegments();
376  if (well->isMultiSegmented()) {
377  start_segment += nseg;
378  start_perforation += nperf;
379  continue;
380  }
381  const Vector cdp_well = subset(cdp, Span(nperf, 1, start_perforation));
382  const ADB segp = subset(state.segp, Span(nseg, 1, start_segment));
383  const Vector perfpressure = (well->wellOps().s2p * segp.value().matrix()).array() + cdp_well;
384  std::copy(perfpressure.data(), perfpressure.data() + nperf, &xw.perfPress()[start_perforation]);
385 
386  start_segment += nseg;
387  start_perforation += nperf;
388  }
389  assert(start_segment == nseg_total_);
390  assert(start_perforation == nperf_total_);
391  }
392 
393 
394 
395 
396 
397  template <class SolutionState>
398  void
399  MultisegmentWells::
400  computeSegmentFluidProperties(const SolutionState& state)
401  {
402  const int np = numPhases();
403  const int nw = msWells().size();
404  const int nseg_total = nseg_total_;
405 
406  if ( !wellOps().has_multisegment_wells ){
407  // not sure if this is needed actually
408  // TODO: to check later if this is really necessary.
409  well_segment_densities_ = ADB::constant(Vector::Zero(nseg_total));
410  segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
411  segment_viscosities_ = ADB::constant(Vector::Zero(nseg_total));
412  for (int phase = 0; phase < np; ++phase) {
413  segment_comp_surf_volume_current_[phase] = ADB::constant(Vector::Zero(nseg_total));
414  segmentCompSurfVolumeInitial()[phase] = Vector::Zero(nseg_total);
415  }
416  return;
417  }
418 
419  // although we will calculate segment density for non-segmented wells at the same time,
420  // while under most of the cases, they will not be used,
421  // since for most of the cases, the density calculation for non-segment wells are
422  // set to be 'SEG' way, which is not a option for multi-segment wells.
423  // When the density calcuation for non-segmented wells are set to 'AVG', then
424  // the density calculation of the mixtures can be the same, while it remains to be verified.
425 
426  // The grid cells associated with segments.
427  // TODO: shoud be computed once and stored in WellState or global Wells structure or class.
428  std::vector<int> segment_cells;
429  segment_cells.reserve(nseg_total);
430  for (int w = 0; w < nw; ++w) {
431  const std::vector<int>& segment_cells_well = msWells()[w]->segmentCells();
432  segment_cells.insert(segment_cells.end(), segment_cells_well.begin(), segment_cells_well.end());
433  }
434  assert(int(segment_cells.size()) == nseg_total);
435 
436  const ADB segment_temp = subset(state.temperature, segment_cells);
437  // using the segment pressure or the average pressure
438  // using the segment pressure first
439  const ADB& segment_press = state.segp;
440 
441  // Compute PVT properties for segments.
442  std::vector<PhasePresence> segment_cond(nseg_total);
443  for (int s = 0; s < nseg_total; ++s) {
444  segment_cond[s] = (*phase_condition_)[segment_cells[s]];
445  }
446  std::vector<ADB> b_seg(np, ADB::null());
447  // Viscosities for different phases
448  std::vector<ADB> mu_seg(np, ADB::null());
449  ADB rsmax_seg = ADB::null();
450  ADB rvmax_seg = ADB::null();
451  const PhaseUsage& pu = fluid_->phaseUsage();
452  if (pu.phase_used[Water]) {
453  b_seg[pu.phase_pos[Water]] = fluid_->bWat(segment_press, segment_temp, segment_cells);
454  mu_seg[pu.phase_pos[Water]] = fluid_->muWat(segment_press, segment_temp, segment_cells);
455  }
456  assert((*active_)[Oil]);
457  const ADB segment_so = subset(state.saturation[pu.phase_pos[Oil]], segment_cells);
458  if (pu.phase_used[Oil]) {
459  const ADB segment_rs = subset(state.rs, segment_cells);
460  b_seg[pu.phase_pos[Oil]] = fluid_->bOil(segment_press, segment_temp, segment_rs,
461  segment_cond, segment_cells);
462  // rsmax_seg = fluidRsSat(segment_press, segment_so, segment_cells);
463  rsmax_seg = fluid_->rsSat(segment_press, segment_so, segment_cells);
464  mu_seg[pu.phase_pos[Oil]] = fluid_->muOil(segment_press, segment_temp, segment_rs,
465  segment_cond, segment_cells);
466  }
467  assert((*active_)[Gas]);
468  if (pu.phase_used[Gas]) {
469  const ADB segment_rv = subset(state.rv, segment_cells);
470  b_seg[pu.phase_pos[Gas]] = fluid_->bGas(segment_press, segment_temp, segment_rv,
471  segment_cond, segment_cells);
472  // rvmax_seg = fluidRvSat(segment_press, segment_so, segment_cells);
473  rvmax_seg = fluid_->rvSat(segment_press, segment_so, segment_cells);
474  mu_seg[pu.phase_pos[Gas]] = fluid_->muGas(segment_press, segment_temp, segment_rv,
475  segment_cond, segment_cells);
476  }
477 
478  // Extract segment flow by phase (segqs) and compute total surface rate.
479  ADB tot_surface_rate = ADB::constant(Vector::Zero(nseg_total));
480  std::vector<ADB> segqs(np, ADB::null());
481  for (int phase = 0; phase < np; ++phase) {
482  segqs[phase] = subset(state.segqs, Span(nseg_total, 1, phase * nseg_total));
483  tot_surface_rate += segqs[phase];
484  }
485 
486  // TODO: later this will be implmented as a global mapping
487  std::vector<std::vector<double>> comp_frac(np, std::vector<double>(nseg_total, 0.0));
488  int start_segment = 0;
489  for (int w = 0; w < nw; ++w) {
490  WellMultiSegmentConstPtr well = msWells()[w];
491  const int nseg = well->numberOfSegments();
492  const std::vector<double>& comp_frac_well = well->compFrac();
493  for (int phase = 0; phase < np; ++phase) {
494  for (int s = 0; s < nseg; ++s) {
495  comp_frac[phase][s + start_segment] = comp_frac_well[phase];
496  }
497  }
498  start_segment += nseg;
499  }
500  assert(start_segment == nseg_total);
501 
502  // Compute mix.
503  // 'mix' contains the component fractions under surface conditions.
504  std::vector<ADB> mix(np, ADB::null());
505  for (int phase = 0; phase < np; ++phase) {
506  // initialize to be the compFrac for each well,
507  // then update only the one with non-zero total volume rate
508  mix[phase] = ADB::constant(Eigen::Map<Vector>(comp_frac[phase].data(), nseg_total));
509  }
510  // There should be a better way to do this.
511  Selector<double> non_zero_tot_rate(tot_surface_rate.value(), Selector<double>::NotEqualZero);
512  for (int phase = 0; phase < np; ++phase) {
513  mix[phase] = non_zero_tot_rate.select(segqs[phase] / tot_surface_rate, mix[phase]);
514  }
515 
516  // Calculate rs and rv.
517  ADB rs = ADB::constant(Vector::Zero(nseg_total));
518  ADB rv = rs;
519  const int gaspos = pu.phase_pos[Gas];
520  const int oilpos = pu.phase_pos[Oil];
521  Selector<double> non_zero_mix_oilpos(mix[oilpos].value(), Selector<double>::GreaterZero);
522  Selector<double> non_zero_mix_gaspos(mix[gaspos].value(), Selector<double>::GreaterZero);
523  // What is the better way to do this?
524  // big values should not be necessary
525  ADB big_values = ADB::constant(Vector::Constant(nseg_total, 1.e100));
526  ADB mix_gas_oil = non_zero_mix_oilpos.select(mix[gaspos] / mix[oilpos], big_values);
527  ADB mix_oil_gas = non_zero_mix_gaspos.select(mix[oilpos] / mix[gaspos], big_values);
528  if ((*active_)[Oil]) {
529  Vector selectorUnderRsmax = Vector::Zero(nseg_total);
530  Vector selectorAboveRsmax = Vector::Zero(nseg_total);
531  for (int s = 0; s < nseg_total; ++s) {
532  if (mix_gas_oil.value()[s] > rsmax_seg.value()[s]) {
533  selectorAboveRsmax[s] = 1.0;
534  } else {
535  selectorUnderRsmax[s] = 1.0;
536  }
537  }
538  rs = non_zero_mix_oilpos.select(selectorAboveRsmax * rsmax_seg + selectorUnderRsmax * mix_gas_oil, rs);
539  }
540  if ((*active_)[Gas]) {
541  Vector selectorUnderRvmax = Vector::Zero(nseg_total);
542  Vector selectorAboveRvmax = Vector::Zero(nseg_total);
543  for (int s = 0; s < nseg_total; ++s) {
544  if (mix_oil_gas.value()[s] > rvmax_seg.value()[s]) {
545  selectorAboveRvmax[s] = 1.0;
546  } else {
547  selectorUnderRvmax[s] = 1.0;
548  }
549  }
550  rv = non_zero_mix_gaspos.select(selectorAboveRvmax * rvmax_seg + selectorUnderRvmax * mix_oil_gas, rv);
551  }
552 
553  // Calculate the phase fraction under reservoir conditions.
554  std::vector<ADB> x(np, ADB::null());
555  for (int phase = 0; phase < np; ++phase) {
556  x[phase] = mix[phase];
557  }
558  if ((*active_)[Gas] && (*active_)[Oil]) {
559  x[gaspos] = (mix[gaspos] - mix[oilpos] * rs) / (Vector::Ones(nseg_total) - rs * rv);
560  x[oilpos] = (mix[oilpos] - mix[gaspos] * rv) / (Vector::Ones(nseg_total) - rs * rv);
561  }
562 
563  // Compute total reservoir volume to surface volume ratio.
564  ADB volrat = ADB::constant(Vector::Zero(nseg_total));
565  for (int phase = 0; phase < np; ++phase) {
566  volrat += x[phase] / b_seg[phase];
567  }
568 
569  // Compute segment densities.
570  ADB dens = ADB::constant(Vector::Zero(nseg_total));
571  for (int phase = 0; phase < np; ++phase) {
572  const Vector surface_density = fluid_->surfaceDensity(phase, segment_cells);
573  dens += surface_density * mix[phase];
574  }
575  well_segment_densities_ = dens / volrat;
576 
577  // Calculating the surface volume of each component in the segment
578  assert(np == int(segment_comp_surf_volume_current_.size()));
579  const ADB segment_surface_volume = segvdt_ / volrat;
580  for (int phase = 0; phase < np; ++phase) {
581  segment_comp_surf_volume_current_[phase] = segment_surface_volume * mix[phase];
582  }
583 
584  // Mass flow rate of the segments
585  segment_mass_flow_rates_ = ADB::constant(Vector::Zero(nseg_total));
586  for (int phase = 0; phase < np; ++phase) {
587  // TODO: how to remove one repeated surfaceDensity()
588  const Vector surface_density = fluid_->surfaceDensity(phase, segment_cells);
589  segment_mass_flow_rates_ += surface_density * segqs[phase];
590  }
591 
592  // Viscosity of the fluid mixture in the segments
593  segment_viscosities_ = ADB::constant(Vector::Zero(nseg_total));
594  for (int phase = 0; phase < np; ++phase) {
595  segment_viscosities_ += x[phase] * mu_seg[phase];
596  }
597  }
598 
599 
600 
601 
602 
603  template <class SolutionState>
604  void
605  MultisegmentWells::
606  addWellFluxEq(const std::vector<ADB>& cq_s,
607  const SolutionState& state,
608  LinearisedBlackoilResidual& residual)
609  {
610  if ( !localWellsActive() ) {
611  return;
612  }
613  // the well flux equations are for each segment and each phase.
614  // /delta m_p_n / dt - /sigma Q_pi - /sigma q_pj + Q_pn = 0
615  // 1. It is the gain of the amount of the component p in the segment n during the
616  // current time step under stock-tank conditions.
617  // It is used to handle the volume storage effects of the wellbore.
618  // We need the information from the previous step and the crrent time step.
619  // 2. for the second term, it is flow into the segment from the inlet segments,
620  // which are unknown and treated implictly.
621  // 3. for the third term, it is the inflow through the perforations.
622  // 4. for the last term, it is the outlet rates and also the segment rates,
623  // which are the primary variable.
624  const int np = numPhases();
625  const int nseg_total = nseg_total_;
626 
627  ADB segqs = state.segqs;
628 
629  std::vector<ADB> segment_volume_change_dt(np, ADB::null());
630  for (int phase = 0; phase < np; ++phase) {
631  if ( wellOps().has_multisegment_wells ) {
632  // Gain of the surface volume of each component in the segment by dt
633  segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] -
634  segmentCompSurfVolumeInitial()[phase];
635 
636  // Special handling for when we are called from solveWellEq().
637  // TODO: restructure to eliminate need for special treatmemt.
638  if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) {
639  assert(segment_volume_change_dt[phase].numBlocks() > 2);
640  assert(segqs.numBlocks() == 2);
641  segment_volume_change_dt[phase] = wellhelpers::onlyWellDerivs(segment_volume_change_dt[phase]);
642  assert(segment_volume_change_dt[phase].numBlocks() == 2);
643  }
644 
645  const ADB cq_s_seg = wellOps().p2s * cq_s[phase];
646  const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total));
647  segqs -= superset(cq_s_seg + wellOps().s2s_inlets * segqs_phase + segment_volume_change_dt[phase],
648  Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
649  } else {
650  segqs -= superset(wellOps().p2s * cq_s[phase], Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
651  }
652  }
653 
654  residual.well_flux_eq = segqs;
655  }
656 
657 
658 
659 
660 
661  template <class SolutionState, class WellState>
662  void
663  MultisegmentWells::
664  addWellControlEq(const SolutionState& state,
665  const WellState& xw,
666  const Vector& aliveWells,
667  LinearisedBlackoilResidual& residual)
668  {
669  // the name of the function is a a little misleading.
670  // Basically it is the function for the pressure equation.
671  // And also, it work as the control equation when it is the segment
672  if( msWells().empty() ) return;
673 
674  const int np = numPhases();
675  const int nw = msWells().size();
676  const int nseg_total = nseg_total_;
677 
678  ADB aqua = ADB::constant(Vector::Zero(nseg_total));
679  ADB liquid = ADB::constant(Vector::Zero(nseg_total));
680  ADB vapour = ADB::constant(Vector::Zero(nseg_total));
681 
682  if ((*active_)[Water]) {
683  aqua += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Aqua * nseg_total));
684  }
685  if ((*active_)[Oil]) {
686  liquid += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Liquid * nseg_total));
687  }
688  if ((*active_)[Gas]) {
689  vapour += subset(state.segqs, Span(nseg_total, 1, BlackoilPhases::Vapour * nseg_total));
690  }
691 
692  // THP control is not implemented for the moment.
693 
694  // Hydrostatic correction variables
695  Vector rho_v = Vector::Zero(nw);
696  Vector vfp_ref_depth_v = Vector::Zero(nw);
697 
698  // Target vars
699  Vector bhp_targets = Vector::Zero(nw);
700  Vector rate_targets = Vector::Zero(nw);
701  Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
702 
703  // Selection variables
704  // well selectors
705  std::vector<int> bhp_well_elems;
706  std::vector<int> rate_well_elems;
707  // segment selectors
708  std::vector<int> bhp_top_elems;
709  std::vector<int> rate_top_elems;
710  std::vector<int> rate_top_phase_elems;
711  std::vector<int> others_elems;
712 
713  //Run through all wells to calculate BHP/RATE targets
714  //and gather info about current control
715  int start_segment = 0;
716  for (int w = 0; w < nw; ++w) {
717  const struct WellControls* wc = msWells()[w]->wellControls();
718 
719  // The current control in the well state overrides
720  // the current control set in the Wells struct, which
721  // is instead treated as a default.
722  const int current = xw.currentControls()[w];
723 
724  const int nseg = msWells()[w]->numberOfSegments();
725 
726  switch (well_controls_iget_type(wc, current)) {
727  case BHP:
728  {
729  bhp_well_elems.push_back(w);
730  bhp_top_elems.push_back(start_segment);
731  bhp_targets(w) = well_controls_iget_target(wc, current);
732  rate_targets(w) = -1e100;
733  for (int p = 0; p < np; ++p) {
734  rate_top_phase_elems.push_back(np * start_segment + p);
735  }
736  }
737  break;
738 
739  case THP:
740  {
741  OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
742  }
743  break;
744 
745  case RESERVOIR_RATE: // Intentional fall-through
746  case SURFACE_RATE:
747  {
748  rate_well_elems.push_back(w);
749  rate_top_elems.push_back(start_segment);
750  for (int p = 0; p < np; ++p) {
751  rate_top_phase_elems.push_back(np * start_segment + p);
752  }
753  // RESERVOIR and SURFACE rates look the same, from a
754  // high-level point of view, in the system of
755  // simultaneous linear equations.
756 
757  const double* const distr =
758  well_controls_iget_distr(wc, current);
759 
760  for (int p = 0; p < np; ++p) {
761  rate_distr.insert(w, p*nw + w) = distr[p];
762  }
763 
764  bhp_targets(w) = -1.0e100;
765  rate_targets(w) = well_controls_iget_target(wc, current);
766  }
767  break;
768 
769  }
770 
771  for (int i = 1; i < nseg; ++i) {
772  others_elems.push_back(i + start_segment);
773  }
774  start_segment += nseg;
775  }
776 
777  // for each segment: 1, if the segment is the top segment, then control equation
778  // 2, if the segment is not the top segment, then the pressure equation
779  const ADB bhp_residual = subset(state.segp, bhp_top_elems) - subset(bhp_targets, bhp_well_elems);
780  const ADB rate_residual = subset(rate_distr * subset(state.segqs, rate_top_phase_elems) - rate_targets, rate_well_elems);
781 
782  ADB others_residual = ADB::constant(Vector::Zero(nseg_total));
783 
784  if ( wellOps().has_multisegment_wells ) {
785  // Special handling for when we are called from solveWellEq().
786  // TODO: restructure to eliminate need for special treatmemt.
787  ADB wspd = (state.segp.numBlocks() == 2)
788  ? wellhelpers::onlyWellDerivs(well_segment_pressures_delta_)
789  : well_segment_pressures_delta_;
790 
791  others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp + wspd);
792  } else {
793  others_residual = wellOps().eliminate_topseg * (state.segp - wellOps().s2s_outlet * state.segp);
794  }
795 
796  // all the control equations
797  // TODO: can be optimized better
798  ADB well_eq_topsegment = subset(superset(bhp_residual, bhp_top_elems, nseg_total) +
799  superset(rate_residual, rate_top_elems, nseg_total), top_well_segments_);
800 
801  // For wells that are dead (not flowing), and therefore not communicating
802  // with the reservoir, we set the equation to be equal to the well's total
803  // flow. This will be a solution only if the target rate is also zero.
804  Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
805  for (int w = 0; w < nw; ++w) {
806  for (int phase = 0; phase < np; ++phase) {
807  rate_summer.insert(w, phase*nw + w) = 1.0;
808  }
809  }
810  Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
811  // TODO: Here only handles the wells, or the top segments
812  // should we also handle some non-alive non-top segments?
813  // should we introduce the cocept of non-alive segments?
814  // At the moment, we only handle the control equations
815  well_eq_topsegment = alive_selector.select(well_eq_topsegment, rate_summer * subset(state.segqs, rate_top_phase_elems));
816 
817  /* residual_.well_eq = superset(bhp_residual, bhp_top_elems, nseg_total) +
818  superset(rate_residual, rate_top_elems, nseg_total) +
819  superset(others_residual, others_elems, nseg_total); */
820  residual.well_eq = superset(well_eq_topsegment, top_well_segments_, nseg_total) +
821  others_residual;
822  }
823 
824 
825 
826 
827 
828  template <class WellState>
829  void
830  MultisegmentWells::
831  updateWellControls(WellState& xw) const
832  {
833  wellhelpers::WellSwitchingLogger logger;
834 
835  if( msWells().empty() ) return ;
836 
837  // Find, for each well, if any constraints are broken. If so,
838  // switch control to first broken constraint.
839  const int np = numPhases();
840  const int nw = msWells().size();
841  for (int w = 0; w < nw; ++w) {
842  const WellControls* wc = msWells()[w]->wellControls();
843  // The current control in the well state overrides
844  // the current control set in the Wells struct, which
845  // is instead treated as a default.
846  int current = xw.currentControls()[w];
847  // Loop over all controls except the current one, and also
848  // skip any RESERVOIR_RATE controls, since we cannot
849  // handle those.
850  const int nwc = well_controls_get_num(wc);
851  int ctrl_index = 0;
852  for (; ctrl_index < nwc; ++ctrl_index) {
853  if (ctrl_index == current) {
854  // This is the currently used control, so it is
855  // used as an equation. So this is not used as an
856  // inequality constraint, and therefore skipped.
857  continue;
858  }
859  if (wellhelpers::constraintBroken(
860  xw.bhp(), xw.thp(), xw.wellRates(),
861  w, np, msWells()[w]->wellType(), wc, ctrl_index)) {
862  // ctrl_index will be the index of the broken constraint after the loop.
863  break;
864  }
865  }
866 
867  if (ctrl_index != nwc) {
868  // Constraint number ctrl_index was broken, switch to it.
869  // Each well is only active on one process. Therefore we always print the sitch info.
870  logger.wellSwitched(msWells()[w]->name(),
871  well_controls_iget_type(wc, current),
872  well_controls_iget_type(wc, ctrl_index));
873  xw.currentControls()[w] = ctrl_index;
874  current = xw.currentControls()[w];
875  }
876 
877  // Get gravity for THP hydrostatic corrrection
878  // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
879 
880  // Updating well state and primary variables.
881  // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
882  const double target = well_controls_iget_target(wc, current);
883  const double* distr = well_controls_iget_distr(wc, current);
884  switch (well_controls_iget_type(wc, current)) {
885  case BHP:
886  xw.bhp()[w] = target;
887  xw.segPress()[top_well_segments_[w]] = target;
888  break;
889 
890  case THP: {
891  OPM_THROW(std::runtime_error, "THP control is not implemented for multi-sgement wells yet!!");
892  }
893 
894  case RESERVOIR_RATE:
895  // No direct change to any observable quantity at
896  // surface condition. In this case, use existing
897  // flow rates as initial conditions as reservoir
898  // rate acts only in aggregate.
899  break;
900 
901  case SURFACE_RATE:
902  for (int phase = 0; phase < np; ++phase) {
903  if (distr[phase] > 0.0) {
904  xw.wellRates()[np * w + phase] = target * distr[phase];
905  // TODO: consider changing all (not just top) segment rates
906  // to make them consistent, it could possibly improve convergence.
907  xw.segPhaseRates()[np * xw.topSegmentLoc()[w] + phase] = target * distr[phase];
908  }
909  }
910  break;
911  }
912 
913  if (wellCollection()->groupControlActive()) {
914  // get the well node in the well collection
915  WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
916 
917  // update whehter the well is under group control or individual control
918  if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
919  // under group control
920  well_node.setIndividualControl(false);
921  } else {
922  // individual control
923  well_node.setIndividualControl(true);
924  }
925  }
926 
927  }
928 
929  }
930 
931 
932 
933 
934 
935  // TODO: This is just a preliminary version, remains to be improved later when we decide a better way
936  // TODO: to intergrate the usual wells and multi-segment wells.
937  template <class SolutionState, class WellState>
938  void
939  MultisegmentWells::
940  computeWellConnectionPressures(const SolutionState& state,
941  const WellState& xw,
942  const std::vector<ADB>& kr_adb,
943  const std::vector<ADB>& fluid_density)
944  {
945  if( ! wellsActive() ) return ;
946 
947  using namespace Opm::AutoDiffGrid;
948  // 1. Compute properties required by computeConnectionPressureDelta().
949  // Note that some of the complexity of this part is due to the function
950  // taking std::vector<double> arguments, and not Eigen objects.
951  const int nperf_total = nperf_total_;
952  const int nw = numWells();
953 
954  const std::vector<int>& well_cells = wellOps().well_cells;
955 
956  well_perforation_densities_ = Vector::Zero(nperf_total);
957 
958  const Vector perf_press = Eigen::Map<const Vector>(xw.perfPress().data(), nperf_total);
959 
960  Vector avg_press = perf_press * 0.0;
961 
962  // for the non-segmented/regular wells, calculated the average pressures.
963  // If it is the top perforation, then average with the bhp().
964  // If it is not the top perforation, then average with the perforation above it().
965  int start_segment = 0;
966  for (int w = 0; w < nw; ++w) {
967  const int nseg = wells_multisegment_[w]->numberOfSegments();
968  if (wells_multisegment_[w]->isMultiSegmented()) {
969  // maybe we should give some reasonable values to prevent the following calculations fail
970  start_segment += nseg;
971  continue;
972  }
973 
974  std::string well_name(wells_multisegment_[w]->name());
975  typedef typename WellState::SegmentedWellMapType::const_iterator const_iterator;
976  const_iterator it_well = xw.segmentedWellMap().find(well_name);
977  assert(it_well != xw.segmentedWellMap().end());
978 
979  const int start_perforation = (*it_well).second.start_perforation;
980  const int end_perforation = start_perforation + (*it_well).second.number_of_perforations;
981  for (int perf = start_perforation; perf < end_perforation; ++perf) {
982  const double p_above = perf == start_perforation ? state.segp.value()[start_segment] : perf_press[perf - 1];
983  const double p_avg = (perf_press[perf] + p_above)/2;
984  avg_press[perf] = p_avg;
985  }
986  start_segment += nseg;
987  }
988  assert(start_segment == xw.numSegments());
989 
990  // Use cell values for the temperature as the wells don't knows its temperature yet.
991  const ADB perf_temp = subset(state.temperature, well_cells);
992 
993  // Compute b, rsmax, rvmax values for perforations.
994  // Evaluate the properties using average well block pressures
995  // and cell values for rs, rv, phase condition and temperature.
996  const ADB avg_press_ad = ADB::constant(avg_press);
997  std::vector<PhasePresence> perf_cond(nperf_total);
998  for (int perf = 0; perf < nperf_total; ++perf) {
999  perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
1000  }
1001  const PhaseUsage& pu = fluid_->phaseUsage();
1002  DataBlock b(nperf_total, pu.num_phases);
1003  std::vector<double> rsmax_perf(nperf_total, 0.0);
1004  std::vector<double> rvmax_perf(nperf_total, 0.0);
1005  if (pu.phase_used[BlackoilPhases::Aqua]) {
1006  const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
1007  b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
1008  }
1009  assert((*active_)[Oil]);
1010  const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
1011  if (pu.phase_used[BlackoilPhases::Liquid]) {
1012  const ADB perf_rs = subset(state.rs, well_cells);
1013  const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
1014  b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
1015  const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
1016  rsmax_perf.assign(rssat.data(), rssat.data() + nperf_total);
1017  }
1018  if (pu.phase_used[BlackoilPhases::Vapour]) {
1019  const ADB perf_rv = subset(state.rv, well_cells);
1020  const Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
1021  b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
1022  const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
1023  rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf_total);
1024  }
1025  // b is row major, so can just copy data.
1026  std::vector<double> b_perf(b.data(), b.data() + nperf_total * pu.num_phases);
1027 
1028  const Vector& perfcelldepth = perf_cell_depth_;
1029  std::vector<double> perf_cell_depth(perfcelldepth.data(), perfcelldepth.data() + nperf_total);
1030 
1031  // Surface density.
1032  // The compute density segment wants the surface densities as
1033  // an np * number of wells cells array
1034  Vector rho = superset(fluid_->surfaceDensity(0 , well_cells), Span(nperf_total, pu.num_phases, 0), nperf_total * pu.num_phases);
1035  for (int phase = 1; phase < pu.num_phases; ++phase) {
1036  rho += superset(fluid_->surfaceDensity(phase , well_cells), Span(nperf_total, pu.num_phases, phase), nperf_total * pu.num_phases);
1037  }
1038  std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf_total * pu.num_phases);
1039 
1040  // 2. Compute densities
1041  std::vector<double> cd =
1043  wells(), fluid_->phaseUsage(), xw.perfPhaseRates(),
1044  b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
1045 
1046  // 3. Compute pressure deltas
1047  std::vector<double> cdp =
1049  wells(), perf_cell_depth, cd, gravity_);
1050 
1051  // 4. Store the results
1052  well_perforation_densities_ = Eigen::Map<const Vector>(cd.data(), nperf_total); // This one is not useful for segmented wells at all
1053  well_perforation_pressure_diffs_ = Eigen::Map<const Vector>(cdp.data(), nperf_total);
1054 
1055  if ( !wellOps().has_multisegment_wells ) {
1056  well_perforation_cell_densities_ = Vector::Zero(nperf_total);
1057  well_perforation_cell_pressure_diffs_ = Vector::Zero(nperf_total);
1058  return;
1059  }
1060 
1061  // compute the average of the fluid densites in the well blocks.
1062  // the average is weighted according to the fluid relative permeabilities.
1063  // const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
1064  size_t temp_size = kr_adb.size();
1065  std::vector<Vector> perf_kr;
1066  for(size_t i = 0; i < temp_size; ++i) {
1067  // const ADB kr_phase_adb = subset(kr_adb[i], well_cells);
1068  const Vector kr_phase = (subset(kr_adb[i], well_cells)).value();
1069  perf_kr.push_back(kr_phase);
1070  }
1071 
1072 
1073  // compute the averaged density for the well block
1074  // TODO: for the non-segmented wells, they should be set to zero
1075  // TODO: for the moment, they are still calculated, while not used later.
1076  for (int i = 0; i < nperf_total; ++i) {
1077  double sum_kr = 0.;
1078  int np = perf_kr.size(); // make sure it is 3
1079  for (int p = 0; p < np; ++p) {
1080  sum_kr += perf_kr[p][i];
1081  }
1082 
1083  for (int p = 0; p < np; ++p) {
1084  perf_kr[p][i] /= sum_kr;
1085  }
1086  }
1087 
1088  Vector rho_avg_perf = Vector::Constant(nperf_total, 0.0);
1089  // TODO: make sure the order of the density and the order of the kr are the same.
1090  for (int phaseIdx = 0; phaseIdx < fluid_->numPhases(); ++phaseIdx) {
1091  // const int canonicalPhaseIdx = canph_[phaseIdx];
1092  // const ADB fluid_density = fluidDensity(canonicalPhaseIdx, rq_[phaseIdx].b, state.rs, state.rv);
1093  const Vector rho_perf = subset(fluid_density[phaseIdx], well_cells).value();
1094  // TODO: phaseIdx or canonicalPhaseIdx ?
1095  rho_avg_perf += rho_perf * perf_kr[phaseIdx];
1096  }
1097 
1098  well_perforation_cell_densities_ = Eigen::Map<const Vector>(rho_avg_perf.data(), nperf_total);
1099 
1100  well_perforation_cell_pressure_diffs_ = gravity_ * well_perforation_cell_densities_ * perf_cell_depth_diffs_;
1101  }
1102 
1103 
1104 
1105 
1106 
1107  template <class SolutionState>
1108  void
1109  MultisegmentWells::
1110  variableStateExtractWellsVars(const std::vector<int>& indices,
1111  std::vector<ADB>& vars,
1112  SolutionState& state) const
1113  {
1114  // TODO: using the original Qs for the segment rates for now, to be fixed eventually.
1115  // TODO: using the original Bhp for the segment pressures for now, to be fixed eventually.
1116 
1117  // segment phase rates in surface volume
1118  state.segqs = std::move(vars[indices[Qs]]);
1119 
1120  // segment pressures
1121  state.segp = std::move(vars[indices[Bhp]]);
1122 
1123  // The qs and bhp are no longer primary variables, but could
1124  // still be used in computations. They are identical to the
1125  // pressures and flows of the top segments.
1126  const int np = num_phases_;
1127  const int ns = nseg_total_;
1128  const int nw = numWells();
1129  state.qs = ADB::constant(Vector::Zero(np * nw));
1130  for (int phase = 0; phase < np; ++phase) {
1131  // Extract segment fluxes for this phase (ns consecutive elements).
1132  ADB segqs_phase = subset(state.segqs, Span(ns, 1, ns*phase));
1133  // Extract top segment fluxes (= well fluxes)
1134  ADB wellqs_phase = subset(segqs_phase, topWellSegments());
1135  // Expand to full size of qs (which contains all phases) and add.
1136  state.qs += superset(wellqs_phase, Span(nw, 1, nw * phase), nw * np);
1137  }
1138  state.bhp = subset(state.segp, topWellSegments());
1139  }
1140 
1141 
1142 
1143 
1144 
1145  template <class WellState>
1146  void
1147  MultisegmentWells::
1148  variableWellStateInitials(const WellState& xw,
1149  std::vector<Vector>& vars0) const
1150  {
1151  // Initial well rates
1152  if ( wells_multisegment_.size() > 0 )
1153  {
1154  // Need to reshuffle well segment rates, from phase running fastest
1155  const int nseg = xw.numSegments();
1156  const int np = xw.numPhases();
1157 
1158  // The transpose() below switches the ordering of the segment rates
1159  const DataBlock segrates = Eigen::Map<const DataBlock>(& xw.segPhaseRates()[0], nseg, np).transpose();
1160  // segment phase rates in surface volume
1161  const Vector segqs = Eigen::Map<const Vector>(segrates.data(), nseg * np);
1162  vars0.push_back(segqs);
1163 
1164  // for the pressure of the segments
1165  const Vector segp = Eigen::Map<const Vector>(& xw.segPress()[0], xw.segPress().size());
1166  vars0.push_back(segp);
1167  }
1168  else
1169  {
1170  // push null sates for segqs and segp
1171  vars0.push_back(Vector());
1172  vars0.push_back(Vector());
1173  }
1174  }
1175 
1176 
1177 }
1178 #endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
ADB muOil(const ADB &po, const ADB &T, const ADB &rs, const std::vector< PhasePresence > &cond, const Cells &cells) const
Oil viscosity.
Definition: BlackoilPropsAdFromDeck.cpp:314
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
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
V surfaceDensity(const int phaseIdx, const Cells &cells) const
Densities of stock components at surface conditions.
Definition: BlackoilPropsAdFromDeck.cpp:242
ADB muWat(const ADB &pw, const ADB &T, const Cells &cells) const
Water viscosity.
Definition: BlackoilPropsAdFromDeck.cpp:263
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
ADB rvSat(const ADB &po, const Cells &cells) const
Condensation curve for Rv as function of oil pressure.
Definition: BlackoilPropsAdFromDeck.cpp:688
V rsSat(const V &po, const Cells &cells) const
Bubble point curve for Rs as function of oil pressure.
Definition: GridHelpers.cpp:27
PhaseUsage phaseUsage() const
Definition: BlackoilPropsAdFromDeck.cpp:231
ADB bWat(const ADB &pw, const ADB &T, const Cells &cells) const
Water formation volume factor.
Definition: BlackoilPropsAdFromDeck.cpp:446
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
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 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
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
ADB muGas(const ADB &pg, const ADB &T, const ADB &rv, const std::vector< PhasePresence > &cond, const Cells &cells) const
Gas viscosity.
Definition: BlackoilPropsAdFromDeck.cpp:381
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