StandardWellsSolvent_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/StandardWellsSolvent.hpp>
23 
24 
25 
26 
27 namespace Opm
28 {
29 
30 
31 
32 
33 
34 
35  StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, WellCollection* well_collection)
36  : Base(wells_arg, well_collection)
37  , solvent_props_(nullptr)
38  , solvent_pos_(-1)
39  , has_solvent_(false)
40  {
41  }
42 
43 
44 
45 
46 
47  void
48  StandardWellsSolvent::initSolvent(const SolventPropsAdFromDeck* solvent_props,
49  const int solvent_pos,
50  const bool has_solvent)
51  {
52  solvent_props_ = solvent_props;
53  solvent_pos_ = solvent_pos;
54  has_solvent_ = has_solvent;
55  }
56 
57 
58 
59 
60 
61  template<class SolutionState, class WellState>
62  void
63  StandardWellsSolvent::
64  computePropertiesForWellConnectionPressures(const SolutionState& state,
65  const WellState& xw,
66  std::vector<double>& b_perf,
67  std::vector<double>& rsmax_perf,
68  std::vector<double>& rvmax_perf,
69  std::vector<double>& surf_dens_perf)
70  {
71  // 1. Compute properties required by computeConnectionPressureDelta().
72  // Note that some of the complexity of this part is due to the function
73  // taking std::vector<double> arguments, and not Eigen objects.
74  const int nperf = wells().well_connpos[wells().number_of_wells];
75  const int nw = wells().number_of_wells;
76 
77  // Compute the average pressure in each well block
78  const Vector perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
79  Vector avg_press = perf_press*0;
80  for (int w = 0; w < nw; ++w) {
81  for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
82  const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
83  const double p_avg = (perf_press[perf] + p_above)/2;
84  avg_press[perf] = p_avg;
85  }
86  }
87 
88  const std::vector<int>& well_cells = wellOps().well_cells;
89 
90  // Use cell values for the temperature as the wells don't knows its temperature yet.
91  const ADB perf_temp = subset(state.temperature, well_cells);
92 
93  // Compute b, rsmax, rvmax values for perforations.
94  // Evaluate the properties using average well block pressures
95  // and cell values for rs, rv, phase condition and temperature.
96  const ADB avg_press_ad = ADB::constant(avg_press);
97  std::vector<PhasePresence> perf_cond(nperf);
98  for (int perf = 0; perf < nperf; ++perf) {
99  perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
100  }
101 
102  const PhaseUsage& pu = fluid_->phaseUsage();
103  DataBlock b(nperf, pu.num_phases);
104 
105  const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
106  if (pu.phase_used[BlackoilPhases::Aqua]) {
107  b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
108  }
109 
110  assert((*active_)[Oil]);
111  assert((*active_)[Gas]);
112  const ADB perf_rv = subset(state.rv, well_cells);
113  const ADB perf_rs = subset(state.rs, well_cells);
114  const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
115  if (pu.phase_used[BlackoilPhases::Liquid]) {
116  const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
117  //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value();
118  b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
119  // const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells);
120  const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
121  rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
122  } else {
123  rsmax_perf.assign(0.0, nperf);
124  }
125  V surf_dens_copy = superset(fluid_->surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
126  for (int phase = 1; phase < pu.num_phases; ++phase) {
127  if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
128  continue; // the gas surface density is added after the solvent is accounted for.
129  }
130  surf_dens_copy += superset(fluid_->surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
131  }
132 
133  if (pu.phase_used[BlackoilPhases::Vapour]) {
134  // Unclear wether the effective or the pure values should be used for the wells
135  // the current usage of unmodified properties values gives best match.
136  //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value();
137  Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
138  Vector rhog = fluid_->surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
139  // to handle solvent related
140  if (has_solvent_) {
141 
142  const Vector bs = solvent_props_->bSolvent(avg_press_ad,well_cells).value();
143  //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
144 
145  // number of cells
146  const int nc = state.pressure.size();
147 
148  const ADB zero = ADB::constant(Vector::Zero(nc));
149  const ADB& ss = state.solvent_saturation;
150  const ADB& sg = ((*active_)[ Gas ]
151  ? state.saturation[ pu.phase_pos[ Gas ] ]
152  : zero);
153 
154  Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
155  Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
156 
157  Vector injectedSolventFraction = Eigen::Map<const Vector>(&xw.solventFraction()[0], nperf);
158 
159  Vector isProducer = Vector::Zero(nperf);
160  Vector ones = Vector::Constant(nperf,1.0);
161  for (int w = 0; w < nw; ++w) {
162  if(wells().type[w] == PRODUCER) {
163  for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
164  isProducer[perf] = 1;
165  }
166  }
167  }
168 
169  F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
170 
171  bg = bg * (ones - F_solvent);
172  bg = bg + F_solvent * bs;
173 
174  const Vector& rhos = solvent_props_->solventSurfaceDensity(well_cells);
175  rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
176  }
177  b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
178  surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
179 
180  // const Vector rvsat = fluidRvSat(avg_press, perf_so, well_cells);
181  const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
182  rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
183  } else {
184  rvmax_perf.assign(0.0, nperf);
185  }
186 
187  // b and surf_dens_perf is row major, so can just copy data.
188  b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
189  surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
190  }
191 
192  template <class SolutionState>
193  void
194  StandardWellsSolvent::
195  computeWellFlux(const SolutionState& state,
196  const std::vector<ADB>& mob_perfcells,
197  const std::vector<ADB>& b_perfcells,
198  Vector& aliveWells,
199  std::vector<ADB>& cq_s) const
200  {
201  if( ! localWellsActive() ) return ;
202 
203  const int np = wells().number_of_phases;
204  const int nw = wells().number_of_wells;
205  const int nperf = wells().well_connpos[nw];
206  Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
207  const std::vector<int>& well_cells = wellOps().well_cells;
208 
209  // pressure diffs computed already (once per step, not changing per iteration)
210  const Vector& cdp = wellPerforationPressureDiffs();
211  // Extract needed quantities for the perforation cells
212  const ADB& p_perfcells = subset(state.pressure, well_cells);
213 
214  // Perforation pressure
215  const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
216 
217  // Pressure drawdown (also used to determine direction of flow)
218  const ADB drawdown = p_perfcells - perfpressure;
219 
220  // Compute vectors with zero and ones that
221  // selects the wanted quantities.
222 
223  // selects injection perforations
224  Vector selectInjectingPerforations = Vector::Zero(nperf);
225  // selects producing perforations
226  Vector selectProducingPerforations = Vector::Zero(nperf);
227  for (int c = 0; c < nperf; ++c){
228  if (drawdown.value()[c] < 0)
229  selectInjectingPerforations[c] = 1;
230  else
231  selectProducingPerforations[c] = 1;
232  }
233 
234  // Handle cross flow
235  const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
236  const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
237  for (int w = 0; w < nw; ++w) {
238  if (!wells().allow_cf[w]) {
239  for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
240  // Crossflow is not allowed; reverse flow is prevented.
241  // At least one of the perforation must be open in order to have a meeningful
242  // equation to solve. For the special case where all perforations have reverse flow,
243  // and the target rate is non-zero all of the perforations are keept open.
244  if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
245  selectProducingPerforations[perf] = 0.0;
246  } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
247  selectInjectingPerforations[perf] = 0.0;
248  }
249  }
250  }
251  }
252 
253  // HANDLE FLOW INTO WELLBORE
254  // compute phase volumetric rates at standard conditions
255  std::vector<ADB> cq_p(np, ADB::null());
256  std::vector<ADB> cq_ps(np, ADB::null());
257  for (int phase = 0; phase < np; ++phase) {
258  cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
259  cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
260  }
261 
262  Vector ones = Vector::Constant(nperf,1.0);
263  ADB F_gas = ADB::constant(ones);
264 
265  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
266  if ((*active_)[Oil] && (*active_)[Gas]) {
267  const int oilpos = pu.phase_pos[Oil];
268  const int gaspos = pu.phase_pos[Gas];
269  const ADB cq_psOil = cq_ps[oilpos];
270  ADB cq_psGas = cq_ps[gaspos];
271  const ADB& rv_perfcells = subset(state.rv, well_cells);
272  const ADB& rs_perfcells = subset(state.rs, well_cells);
273  cq_ps[gaspos] += rs_perfcells * cq_psOil;
274 
275  if(has_solvent_) {
276  // The solvent gas need to be removed from the gas
277  // before multiplied with rv.
278  const ADB& ss = state.solvent_saturation;
279  const ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
280 
281  Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
282  F_gas -= subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
283  cq_psGas = cq_psGas * F_gas;
284  }
285  cq_ps[oilpos] += rv_perfcells * cq_psGas;
286  }
287 
288  // HANDLE FLOW OUT FROM WELLBORE
289  // Using total mobilities
290  ADB total_mob = mob_perfcells[0];
291  for (int phase = 1; phase < np; ++phase) {
292  total_mob += mob_perfcells[phase];
293  }
294  // injection perforations total volume rates
295  const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
296 
297  // Store well perforation total fluxes (reservor volumes) if requested.
298  if (store_well_perforation_fluxes_) {
299  // Ugly const-cast, but unappealing alternatives.
300  Vector& wf = const_cast<Vector&>(well_perforation_fluxes_);
301  wf = cqt_i.value();
302  for (int phase = 0; phase < np; ++phase) {
303  wf += cq_p[phase].value();
304  }
305  }
306 
307  // compute wellbore mixture for injecting perforations
308  // The wellbore mixture depends on the inflow from the reservoar
309  // and the well injection rates.
310 
311  // compute avg. and total wellbore phase volumetric rates at standard conds
312  const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
313  std::vector<ADB> wbq(np, ADB::null());
314  ADB wbqt = ADB::constant(Vector::Zero(nw));
315  for (int phase = 0; phase < np; ++phase) {
316  const ADB& q_ps = wellOps().p2w * cq_ps[phase];
317  const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
318  Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
319  const int pos = pu.phase_pos[phase];
320  wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps;
321  wbqt += wbq[phase];
322  }
323  // compute wellbore mixture at standard conditions.
324  Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
325  std::vector<ADB> cmix_s(np, ADB::null());
326  for (int phase = 0; phase < np; ++phase) {
327  const int pos = pu.phase_pos[phase];
328  cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
329  }
330 
331  // compute volume ratio between connection at standard conditions
332  ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
333 
334  if ((*active_)[Water]) {
335  const int watpos = pu.phase_pos[Water];
336  volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
337  }
338 
339  if ((*active_)[Oil] && (*active_)[Gas]) {
340  // Incorporate RS/RV factors if both oil and gas active
341  const ADB& rv_perfcells = subset(state.rv, well_cells);
342  const ADB& rs_perfcells = subset(state.rs, well_cells);
343  const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
344 
345  const int oilpos = pu.phase_pos[Oil];
346  const int gaspos = pu.phase_pos[Gas];
347 
348  const ADB tmp_oil = (cmix_s[oilpos] - rv_perfcells * F_gas * cmix_s[gaspos]) / d;
349  volumeRatio += tmp_oil / b_perfcells[oilpos];
350 
351  const ADB tmp_gas = (cmix_s[gaspos] - rs_perfcells * cmix_s[oilpos]) / d;
352  volumeRatio += tmp_gas / b_perfcells[gaspos];
353  }
354  else {
355  if ((*active_)[Oil]) {
356  const int oilpos = pu.phase_pos[Oil];
357  volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
358  }
359  if ((*active_)[Gas]) {
360  const int gaspos = pu.phase_pos[Gas];
361  volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
362  }
363  }
364 
365 
366  // injecting connections total volumerates at standard conditions
367  ADB cqt_is = cqt_i/volumeRatio;
368 
369  // connection phase volumerates at standard conditions
370  cq_s.resize(np, ADB::null());
371  for (int phase = 0; phase < np; ++phase) {
372  cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
373  }
374 
375  // check for dead wells (used in the well controll equations)
376  aliveWells = Vector::Constant(nw, 1.0);
377  for (int w = 0; w < nw; ++w) {
378  if (wbqt.value()[w] == 0) {
379  aliveWells[w] = 0.0;
380  }
381  }
382  }
383 
384 
385 
386 
387 
388 
389  template <class SolutionState, class WellState>
390  void
391  StandardWellsSolvent::
392  computeWellConnectionPressures(const SolutionState& state,
393  const WellState& xw)
394  {
395  if( ! localWellsActive() ) return ;
396  // 1. Compute properties required by computeConnectionPressureDelta().
397  // Note that some of the complexity of this part is due to the function
398  // taking std::vector<double> arguments, and not Eigen objects.
399  std::vector<double> b_perf;
400  std::vector<double> rsmax_perf;
401  std::vector<double> rvmax_perf;
402  std::vector<double> surf_dens_perf;
403  computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
404 
405  const Vector pdepth = perf_cell_depth_;
406  const int nperf = wells().well_connpos[wells().number_of_wells];
407  const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
408 
409  computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
410 
411  }
412 
413 
414 
415 
416 
417  template <class ReservoirResidualQuant, class SolutionState>
418  void
419  StandardWellsSolvent::
420  extractWellPerfProperties(const SolutionState& state,
421  const std::vector<ReservoirResidualQuant>& rq,
422  std::vector<ADB>& mob_perfcells,
423  std::vector<ADB>& b_perfcells) const
424  {
425  Base::extractWellPerfProperties(state, rq, mob_perfcells, b_perfcells);
426  // handle the solvent related
427  if (has_solvent_) {
428  const Opm::PhaseUsage& pu = fluid_->phaseUsage();
429  int gas_pos = pu.phase_pos[Gas];
430  const std::vector<int>& well_cells = wellOps().well_cells;
431  const int nperf = well_cells.size();
432  // Gas and solvent is combinded and solved together
433  // The input in the well equation is then the
434  // total gas phase = hydro carbon gas + solvent gas
435 
436  // The total mobility is the sum of the solvent and gas mobiliy
437  mob_perfcells[gas_pos] += subset(rq[solvent_pos_].mob, well_cells);
438 
439  // A weighted sum of the b-factors of gas and solvent are used.
440  const int nc = rq[solvent_pos_].mob.size();
441 
442  const ADB zero = ADB::constant(Vector::Zero(nc));
443  const ADB& ss = state.solvent_saturation;
444  const ADB& sg = ((*active_)[ Gas ]
445  ? state.saturation[ pu.phase_pos[ Gas ] ]
446  : zero);
447 
448  Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
449  ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
450  Vector ones = Vector::Constant(nperf,1.0);
451 
452  b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
453  b_perfcells[gas_pos] += (F_solvent * subset(rq[solvent_pos_].b, well_cells));
454  }
455  }
456 
457 
458 }
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
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
ADB bSolvent(const ADB &pg, const Cells &cells) const
Solvent formation volume factor.
Definition: SolventPropsAdFromDeck.cpp:339
Vector & wellPerforationPressureDiffs()
Diff to bhp for each well perforation.
Definition: StandardWells_impl.hpp:204
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 solventSurfaceDensity(const Cells &cells) const
Solvent surface density.
Definition: SolventPropsAdFromDeck.cpp:447
V rsSat(const V &po, const Cells &cells) const
Bubble point curve for Rs as function of oil pressure.
bool localWellsActive() const
return true if wells are available on this process
Definition: StandardWells_impl.hpp:148
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
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