BlackoilSolventModel_impl.hpp
1 /*
2  Copyright 2015 IRIS AS
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED
21 #define OPM_BLACKOILSOLVENTMODEL_IMPL_HEADER_INCLUDED
22 
23 #include <opm/autodiff/BlackoilSolventModel.hpp>
24 
25 #include <opm/autodiff/AutoDiffBlock.hpp>
26 #include <opm/autodiff/AutoDiffHelpers.hpp>
27 #include <opm/autodiff/GridHelpers.hpp>
28 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
29 #include <opm/autodiff/GeoProps.hpp>
30 #include <opm/autodiff/WellDensitySegmented.hpp>
31 
32 #include <opm/core/grid.h>
33 #include <opm/core/linalg/LinearSolverInterface.hpp>
34 #include <opm/core/linalg/ParallelIstlInformation.hpp>
35 #include <opm/core/props/rock/RockCompressibility.hpp>
36 #include <opm/common/ErrorMacros.hpp>
37 #include <opm/common/Exceptions.hpp>
38 #include <opm/parser/eclipse/Units/Units.hpp>
39 #include <opm/core/well_controls.h>
40 #include <opm/core/utility/parameters/ParameterGroup.hpp>
41 
42 #include <cassert>
43 #include <cmath>
44 #include <iostream>
45 #include <iomanip>
46 #include <limits>
47 
48 namespace Opm {
49 
50 
51 
52  namespace detail {
53 
54  template <class PU>
55  int solventPos(const PU& pu)
56  {
57  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
58  int pos = 0;
59  for (int phase = 0; phase < maxnp; ++phase) {
60  if (pu.phase_used[phase]) {
61  pos++;
62  }
63  }
64 
65  return pos;
66  }
67 
68  } // namespace detail
69 
70 
71 
72  template <class Grid>
73  BlackoilSolventModel<Grid>::BlackoilSolventModel(const typename Base::ModelParameters& param,
74  const Grid& grid,
75  const BlackoilPropsAdFromDeck& fluid,
76  const DerivedGeology& geo,
77  const RockCompressibility* rock_comp_props,
78  const SolventPropsAdFromDeck& solvent_props,
79  const StandardWellsSolvent& well_model,
80  const NewtonIterationBlackoilInterface& linsolver,
81  std::shared_ptr< const EclipseState > eclState,
82  const bool has_disgas,
83  const bool has_vapoil,
84  const bool terminal_output,
85  const bool has_solvent,
86  const bool is_miscible)
87  : Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver,
88  eclState, has_disgas, has_vapoil, terminal_output),
89  has_solvent_(has_solvent),
90  solvent_pos_(detail::solventPos(fluid.phaseUsage())),
91  solvent_props_(solvent_props),
92  is_miscible_(is_miscible)
93 
94  {
95  if (has_solvent_) {
96 
97  // If deck has solvent, residual_ should contain solvent equation.
98  sd_.rq.resize(fluid_.numPhases() + 1);
99  residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
100  Base::material_name_.push_back("Solvent");
101  assert(solvent_pos_ == fluid_.numPhases());
102  residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
103 
104  wellModel().initSolvent(&solvent_props_, solvent_pos_, has_solvent_);
105  }
106  if (is_miscible_) {
107  mu_eff_.resize(fluid_.numPhases() + 1, ADB::null());
108  b_eff_.resize(fluid_.numPhases() + 1, ADB::null());
109  }
110  }
111 
112 
113 
114 
115 
116  template <class Grid>
117  void
118  BlackoilSolventModel<Grid>::makeConstantState(SolutionState& state) const
119  {
120  Base::makeConstantState(state);
121  state.solvent_saturation = ADB::constant(state.solvent_saturation.value());
122  }
123 
124 
125 
126 
127 
128  template <class Grid>
129  std::vector<V>
130  BlackoilSolventModel<Grid>::variableStateInitials(const ReservoirState& x,
131  const WellState& xw) const
132  {
133  std::vector<V> vars0 = Base::variableStateInitials(x, xw);
134  assert(int(vars0.size()) == fluid_.numPhases() + 2);
135 
136  // Initial solvent concentration.
137  if (has_solvent_) {
138  const auto& solvent_saturation = x.getCellData( BlackoilState::SSOL );
139  const int nc = solvent_saturation.size();
140  const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);
141 
142  // This is some insanely detailed flickety flackety code;
143  // Solvent belongs after other reservoir vars but before well vars.
144  auto solvent_pos = vars0.begin() + fluid_.numPhases();
145  assert (not solvent_saturation.empty());
146  assert(solvent_pos == vars0.end() - 2);
147  vars0.insert(solvent_pos, ss);
148  }
149  return vars0;
150  }
151 
152 
153 
154 
155 
156  template <class Grid>
157  std::vector<int>
158  BlackoilSolventModel<Grid>::variableStateIndices() const
159  {
160  std::vector<int> ind = Base::variableStateIndices();
161  assert(ind.size() == 5);
162  if (has_solvent_) {
163  ind.resize(6);
164  // Solvent belongs after other reservoir vars but before well vars.
165  ind[Solvent] = fluid_.numPhases();
166  // Solvent is pushing back the well vars.
167  ++ind[Qs];
168  ++ind[Bhp];
169  }
170  return ind;
171  }
172 
173 
174 
175 
176  template <class Grid>
177  typename BlackoilSolventModel<Grid>::SolutionState
178  BlackoilSolventModel<Grid>::variableStateExtractVars(const ReservoirState& x,
179  const std::vector<int>& indices,
180  std::vector<ADB>& vars) const
181  {
182  // This is more or less a copy of the base class. Refactoring is needed in the base class
183  // to avoid unnecessary copying.
184 
185  //using namespace Opm::AutoDiffGrid;
186  const int nc = Opm::AutoDiffGrid::numCells(grid_);
187  const Opm::PhaseUsage pu = fluid_.phaseUsage();
188 
189  SolutionState state(fluid_.numPhases());
190 
191  // Pressure.
192  state.pressure = std::move(vars[indices[Pressure]]);
193 
194  // Temperature cannot be a variable at this time (only constant).
195  const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
196  state.temperature = ADB::constant(temp);
197 
198  // Saturations
199  {
200  ADB so = ADB::constant(V::Ones(nc, 1));
201 
202  if (active_[ Water ]) {
203  state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
204  const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
205  so -= sw;
206  }
207  if (has_solvent_) {
208  state.solvent_saturation = std::move(vars[indices[Solvent]]);
209  so -= state.solvent_saturation;
210  }
211 
212  if (active_[ Gas ]) {
213  // Define Sg Rs and Rv in terms of xvar.
214  // Xvar is only defined if gas phase is active
215  const ADB& xvar = vars[indices[Xvar]];
216  ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
217  sg = Base::isSg_*xvar + Base::isRv_*so;
218  so -= sg;
219 
220  if (active_[ Oil ]) {
221  // RS and RV is only defined if both oil and gas phase are active.
222  state.canonical_phase_pressures = computePressures(state.pressure, state.saturation[pu.phase_pos[ Water ]], so, sg, state.solvent_saturation);
223  sd_.rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
224  if (has_disgas_) {
225  state.rs = (1-Base::isRs_)*sd_.rsSat + Base::isRs_*xvar;
226  } else {
227  state.rs = sd_.rsSat;
228  }
229  sd_.rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
230  if (has_vapoil_) {
231  state.rv = (1-Base::isRv_)*sd_.rvSat + Base::isRv_*xvar;
232  } else {
233  state.rv = sd_.rvSat;
234  }
235  }
236  }
237 
238  if (active_[ Oil ]) {
239  // Note that so is never a primary variable.
240  state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
241  }
242  }
243  // wells
244  wellModel().variableStateExtractWellsVars(indices, vars, state);
245  return state;
246 
247  }
248 
249 
250 
251 
252 
253  template <class Grid>
254  void
255  BlackoilSolventModel<Grid>::computeAccum(const SolutionState& state,
256  const int aix )
257  {
258 
259  if (is_miscible_) {
260  computeEffectiveProperties(state);
261  }
262  Base::computeAccum(state, aix);
263 
264  // Compute accumulation of the solvent
265  if (has_solvent_) {
266  const ADB& press = state.pressure;
267  const ADB& ss = state.solvent_saturation;
268  const ADB pv_mult = poroMult(press); // also computed in Base::computeAccum, could be optimized.
269  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
270 
271  const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
272  const std::vector<PhasePresence>& cond = phaseCondition();
273  sd_.rq[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond);
274  sd_.rq[solvent_pos_].accum[aix] = pv_mult * sd_.rq[solvent_pos_].b * ss;
275  }
276  }
277 
278 
279 
280 
281 
282  template <class Grid>
283  void
284  BlackoilSolventModel<Grid>::
285  assembleMassBalanceEq(const SolutionState& state)
286  {
287 
288  Base::assembleMassBalanceEq(state);
289 
290  if (has_solvent_) {
291  residual_.material_balance_eq[ solvent_pos_ ] =
292  pvdt_ * (sd_.rq[solvent_pos_].accum[1] - sd_.rq[solvent_pos_].accum[0])
293  + ops_.div*sd_.rq[solvent_pos_].mflux;
294  }
295 
296  }
297 
298  template <class Grid>
299  void
300  BlackoilSolventModel<Grid>::updateEquationsScaling()
301  {
302  Base::updateEquationsScaling();
303  assert(MaxNumPhases + 1 == residual_.matbalscale.size());
304  if (has_solvent_) {
305  const ADB& temp_b = sd_.rq[solvent_pos_].b;
306  ADB::V B = 1. / temp_b.value();
307 #if HAVE_MPI
308  if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
309  {
310  const ParallelISTLInformation& real_info =
311  boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation());
312  double B_global_sum = 0;
313  real_info.computeReduction(B, Reduction::makeGlobalSumFunctor<double>(), B_global_sum);
314  residual_.matbalscale[solvent_pos_] = B_global_sum / Base::global_nc_;
315  }
316  else
317 #endif
318  {
319  residual_.matbalscale[solvent_pos_] = B.mean();
320  }
321  }
322  }
323 
324  template <class Grid>
325  void BlackoilSolventModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
326  const SolutionState& state,
327  WellState& xw)
328 
329  {
330 
331  // Add well contributions to solvent mass balance equation
332 
333  Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
334 
335  if (has_solvent_) {
336  const int nperf = wells().well_connpos[wells().number_of_wells];
337  const int nc = Opm::AutoDiffGrid::numCells(grid_);
338 
339  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
340  const ADB zero = ADB::constant(V::Zero(nc));
341  const ADB& ss = state.solvent_saturation;
342  const ADB& sg = (active_[ Gas ]
343  ? state.saturation[ pu.phase_pos[ Gas ] ]
344  : zero);
345 
346  const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
347  Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
348  ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
349 
350  const int nw = wells().number_of_wells;
351  V injectedSolventFraction = Eigen::Map<const V>(&xw.solventFraction()[0], nperf);
352 
353  V isProducer = V::Zero(nperf);
354  V ones = V::Constant(nperf,1.0);
355  for (int w = 0; w < nw; ++w) {
356  if(wells().type[w] == PRODUCER) {
357  for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
358  isProducer[perf] = 1;
359  }
360  }
361  }
362 
363  const ADB& rs_perfcells = subset(state.rs, well_cells);
364  const ADB& rv_perfcells = subset(state.rv, well_cells);
365  int gas_pos = fluid_.phaseUsage().phase_pos[Gas];
366  int oil_pos = fluid_.phaseUsage().phase_pos[Oil];
367 
368  // remove contribution from the dissolved gas.
369  const ADB cq_s_solvent = (isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction) * (cq_s[gas_pos] - rs_perfcells * cq_s[oil_pos]) / (ones - rs_perfcells * rv_perfcells);
370 
371  // Solvent contribution to the mass balance equation is given as a fraction
372  // of the gas contribution.
373  residual_.material_balance_eq[solvent_pos_] -= superset(cq_s_solvent, well_cells, nc);
374 
375  // The gas contribution must be reduced accordingly for the total contribution to be
376  // the same.
377  residual_.material_balance_eq[gas_pos] += superset(cq_s_solvent, well_cells, nc);
378 
379  }
380  }
381 
382 
383 
384 
385 
386  template <class Grid>
387  void
389  updateState(const V& dx,
390  ReservoirState& reservoir_state,
391  WellState& well_state)
392  {
393  //TODO:
394  // This is basicly a copy of updateState in the Base class
395  // The convergence is very sensitive to details in the appelyard process
396  // and the hydrocarbonstate detection. Further changes may occur, refactoring
397  // to reuse more of the base class is planned when the code mature a bit more.
398  using namespace Opm::AutoDiffGrid;
399  const int np = fluid_.numPhases();
400  const int nc = numCells(grid_);
401  const V null;
402  assert(null.size() == 0);
403  const V zero = V::Zero(nc);
404  const V ones = V::Constant(nc,1.0);
405 
406  // Extract parts of dx corresponding to each part.
407  const V dp = subset(dx, Span(nc));
408  int varstart = nc;
409  const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null;
410  varstart += dsw.size();
411 
412  const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null;
413  varstart += dxvar.size();
414 
415  const V dss = has_solvent_ ? subset(dx, Span(nc, 1, varstart)) : null;
416  varstart += dss.size();
417 
418  // Extract well parts np phase rates + bhp
419  const V dwells = subset(dx, Span(wellModel().numWellVars(), 1, varstart));
420  varstart += dwells.size();
421 
422  assert(varstart == dx.size());
423 
424  // Pressure update.
425  const double dpmaxrel = dpMaxRel();
426  const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
427  const V absdpmax = dpmaxrel*p_old.abs();
428  const V dp_limited = sign(dp) * dp.abs().min(absdpmax);
429  const V p = (p_old - dp_limited).max(zero);
430  std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
431 
432 
433  // Saturation updates.
434  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
435  const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
436  const double dsmax = dsMax();
437 
438  // initialize with zeros
439  // if the phase is active the saturation are overwritten.
440  V so;
441  V sw = zero;
442  V sg = zero;
443  V ss = zero;
444 
445  // Appleyard chop process.
446  // We chop too large updates in the saturations
447  {
448  V maxVal = zero;
449  V dso = zero;
450  if (active_[Water]){
451  maxVal = dsw.abs().max(maxVal);
452  dso = dso - dsw;
453  }
454 
455  V dsg;
456  if (active_[Gas]){
457  dsg = Base::isSg_ * dxvar - Base::isRv_ * dsw;
458  maxVal = dsg.abs().max(maxVal);
459  dso = dso - dsg;
460  }
461  if (has_solvent_){
462  maxVal = dss.abs().max(maxVal);
463  // solvent is not added note that the so calculated
464  // here is overwritten later
465  //dso = dso - dss;
466  }
467 
468  maxVal = dso.abs().max(maxVal);
469 
470  V step = dsmax/maxVal;
471  step = step.min(1.);
472 
473  if (active_[Water]) {
474  const int pos = pu.phase_pos[ Water ];
475  const V sw_old = s_old.col(pos);
476  sw = sw_old - step * dsw;
477  }
478 
479  if (active_[Gas]) {
480  const int pos = pu.phase_pos[ Gas ];
481  const V sg_old = s_old.col(pos);
482  sg = sg_old - step * dsg;
483  }
484  if (has_solvent_) {
485  auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
486  const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
487  ss = ss_old - step * dss;
488  }
489 
490  const int pos = pu.phase_pos[ Oil ];
491  const V so_old = s_old.col(pos);
492  so = so_old - step * dso;
493  }
494 
495  // solvent is not included in the adjustment for negative saturation
496  auto ixg = sg < 0;
497  for (int c = 0; c < nc; ++c) {
498  if (ixg[c]) {
499  sw[c] = sw[c] / (1-sg[c]);
500  so[c] = so[c] / (1-sg[c]);
501  sg[c] = 0;
502  }
503  }
504 
505  auto ixo = so < 0;
506  for (int c = 0; c < nc; ++c) {
507  if (ixo[c]) {
508  sw[c] = sw[c] / (1-so[c]);
509  sg[c] = sg[c] / (1-so[c]);
510  so[c] = 0;
511  }
512  }
513 
514  auto ixw = sw < 0;
515  for (int c = 0; c < nc; ++c) {
516  if (ixw[c]) {
517  so[c] = so[c] / (1-sw[c]);
518  sg[c] = sg[c] / (1-sw[c]);
519  sw[c] = 0;
520  }
521  }
522 
523  auto ixs = ss < 0;
524  for (int c = 0; c < nc; ++c) {
525  if (ixs[c]) {
526  ss[c] = 0;
527  }
528  }
529 
530 
531  // The oil saturation is defined to
532  // fill the rest of the pore space.
533  // For convergence reasons oil saturations
534  // is included in the appelyard chopping
535  so = V::Constant(nc,1.0) - sw - sg - ss;
536 
537  // Update rs and rv
538  const double drmaxrel = drMaxRel();
539  V rs;
540  if (has_disgas_) {
541  const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
542  const V drs = Base::isRs_ * dxvar;
543  const V drs_limited = sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
544  rs = rs_old - drs_limited;
545  rs = rs.max(zero);
546  }
547  V rv;
548  if (has_vapoil_) {
549  const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
550  const V drv = Base::isRv_ * dxvar;
551  const V drv_limited = sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
552  rv = rv_old - drv_limited;
553  rv = rv.max(zero);
554  }
555 
556  sd_.soMax = fluid_.satOilMax();
557 
558  // Sg is used as primal variable for water only cells.
559  const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
560  auto watOnly = sw > (1 - epsilon);
561 
562  // phase translation sg <-> rs
563  std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
564  std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
565 
566  if (has_disgas_) {
567  const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
568  const V rsSat = fluidRsSat(p, so, cells_);
569  sd_.rsSat = ADB::constant(rsSat);
570  // The obvious case
571  auto hasGas = (sg > 0 && Base::isRs_ == 0);
572 
573  // Set oil saturated if previous rs is sufficiently large
574  const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
575  auto gasVaporized = ( (rs > rsSat * (1+epsilon) && Base::isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
576  auto useSg = watOnly || hasGas || gasVaporized;
577  for (int c = 0; c < nc; ++c) {
578  if (useSg[c]) {
579  rs[c] = rsSat[c];
580  if (watOnly[c]) {
581  so[c] = 0;
582  sg[c] = 0;
583  ss[c] = 0;
584  rs[c] = 0;
585  }
586 
587  } else {
588  hydroCarbonState[c] = HydroCarbonState::OilOnly;
589  }
590  }
591  rs = rs.min(rsSat);
592  }
593 
594  // phase transitions so <-> rv
595  if (has_vapoil_) {
596 
597  // The gas pressure is needed for the rvSat calculations
598  const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
599  const V gaspress = computeGasPressure(p, sw, so, sg);
600  const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
601  const V rvSat = fluidRvSat(gaspress, so, cells_);
602  sd_.rvSat = ADB::constant(rvSat);
603 
604  // The obvious case
605  auto hasOil = (so > 0 && Base::isRv_ == 0);
606 
607  // Set oil saturated if previous rv is sufficiently large
608  const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
609  auto oilCondensed = ( (rv > rvSat * (1+epsilon) && Base::isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
610  auto useSg = watOnly || hasOil || oilCondensed;
611  for (int c = 0; c < nc; ++c) {
612  if (useSg[c]) {
613  rv[c] = rvSat[c];
614  if (watOnly[c]) {
615  so[c] = 0;
616  sg[c] = 0;
617  ss[c] = 0;
618  rv[c] = 0;
619  }
620 
621  } else {
622  hydroCarbonState[c] = HydroCarbonState::GasOnly;
623  }
624  }
625  rv = rv.min(rvSat);
626 
627  }
628 
629  // Update the reservoir_state
630  if (has_solvent_) {
631  auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
632  std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
633  }
634 
635  for (int c = 0; c < nc; ++c) {
636  reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
637  }
638 
639  for (int c = 0; c < nc; ++c) {
640  reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
641  }
642 
643  if (active_[ Oil ]) {
644  const int pos = pu.phase_pos[ Oil ];
645  for (int c = 0; c < nc; ++c) {
646  reservoir_state.saturation()[c*np + pos] = so[c];
647  }
648  }
649 
650  // Update the reservoir_state
651  if (has_disgas_) {
652  std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
653  }
654 
655  if (has_vapoil_) {
656  std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
657  }
658 
659  wellModel().updateWellState(dwells, Base::dbhpMaxRel(), well_state);
660 
661  for( auto w = 0; w < wells().number_of_wells; ++w ) {
662  if (wells().type[w] == INJECTOR) {
663  continue;
664  }
665 
666  for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf ) {
667  int wc = wells().well_cells[perf];
668  if ( (ss[wc] + sg[wc]) > 0) {
669  well_state.solventFraction()[perf] = ss[wc] / (ss[wc] + sg[wc]);
670  } else {
671  well_state.solventFraction()[perf] = 0.0;
672  }
673  }
674  }
675 
676 
677 
678  // Update phase conditions used for property calculations.
679  updatePhaseCondFromPrimalVariable(reservoir_state);
680  }
681 
682 
683 
684  template <class Grid>
685  void
687  const V& transi,
688  const ADB& kr ,
689  const ADB& mu ,
690  const ADB& rho ,
691  const ADB& phasePressure,
692  const SolutionState& state)
693  {
694 
695  const int canonicalPhaseIdx = canph_[ actph ];
696  // make a copy to make it possible to modify it
697  ADB kr_mod = kr;
698  if (canonicalPhaseIdx == Gas) {
699  if (has_solvent_) {
700  const int nc = Opm::UgGridHelpers::numCells(grid_);
701  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
702  const ADB zero = ADB::constant(V::Zero(nc));
703  const V ones = V::Constant(nc, 1.0);
704 
705  const ADB& ss = state.solvent_saturation;
706  const ADB& sg = (active_[ Gas ]
707  ? state.saturation[ pu.phase_pos[ Gas ] ]
708  : zero);
709  Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
710  const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg));
711 
712  // Compute solvent properties
713  const std::vector<PhasePresence>& cond = phaseCondition();
714  ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond);
715  ADB rho_s = fluidDensity(Solvent,sd_.rq[solvent_pos_].b, state.rs, state.rv);
716 
717  // Compute solvent relperm and mass flux
718  ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod;
719  Base::computeMassFlux(solvent_pos_, transi, krs, mu_s, rho_s, phasePressure, state);
720 
721  // Modify gas relperm
722  kr_mod = solvent_props_.gasRelPermMultiplier( (ones - F_solvent) , cells_) * kr_mod;
723 
724  }
725  }
726  // Compute mobility and flux
727  Base::computeMassFlux(actph, transi, kr_mod, mu, rho, phasePressure, state);
728 
729  }
730 
731  template <class Grid>
732  ADB
733  BlackoilSolventModel<Grid>::fluidViscosity(const int phase,
734  const ADB& p ,
735  const ADB& temp ,
736  const ADB& rs ,
737  const ADB& rv ,
738  const std::vector<PhasePresence>& cond) const
739  {
740  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
741  if (phase == Solvent) {
742  if (!is_miscible_) {
743  return solvent_props_.muSolvent(p, cells_);
744  } else {
745  return mu_eff_[solvent_pos_];
746  }
747 
748  } else {
749  if (!is_miscible_) {
750  return Base::fluidViscosity(phase, p, temp, rs, rv, cond);
751  } else {
752  return mu_eff_[pu.phase_pos[ phase ]];
753  }
754  }
755  }
756 
757  template <class Grid>
758  ADB
759  BlackoilSolventModel<Grid>::fluidReciprocFVF(const int phase,
760  const ADB& p ,
761  const ADB& temp ,
762  const ADB& rs ,
763  const ADB& rv ,
764  const std::vector<PhasePresence>& cond) const
765  {
766  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
767  if (phase == Solvent) {
768  if (!is_miscible_) {
769  return solvent_props_.bSolvent(p, cells_);
770  } else {
771  return b_eff_[solvent_pos_];
772  }
773 
774  } else {
775  if (!is_miscible_) {
776  return Base::fluidReciprocFVF(phase, p, temp, rs, rv, cond);
777  } else {
778  return b_eff_[pu.phase_pos[ phase ]];
779  }
780  }
781  }
782 
783  template <class Grid>
784  ADB
785  BlackoilSolventModel<Grid>::fluidDensity(const int phase,
786  const ADB& b,
787  const ADB& rs,
788  const ADB& rv) const
789  {
790  if (phase == Solvent && has_solvent_) {
791  return solvent_props_.solventSurfaceDensity(cells_) * b;
792  } else {
793  return Base::fluidDensity(phase, b, rs, rv);
794  }
795  }
796 
797  template <class Grid>
798  std::vector<ADB>
799  BlackoilSolventModel<Grid>::computeRelPerm(const SolutionState& state) const
800  {
801  using namespace Opm::AutoDiffGrid;
802  const int nc = numCells(grid_);
803 
804  const ADB zero = ADB::constant(V::Zero(nc));
805 
806  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
807  const ADB& sw = (active_[ Water ]
808  ? state.saturation[ pu.phase_pos[ Water ] ]
809  : zero);
810 
811  const ADB& so = (active_[ Oil ]
812  ? state.saturation[ pu.phase_pos[ Oil ] ]
813  : zero);
814 
815  const ADB& sg = (active_[ Gas ]
816  ? state.saturation[ pu.phase_pos[ Gas ] ]
817  : zero);
818 
819  if (has_solvent_) {
820  const ADB& ss = state.solvent_saturation;
821  if (is_miscible_) {
822 
823  assert(active_[ Oil ]);
824  assert(active_[ Gas ]);
825 
826  std::vector<ADB> relperm = fluid_.relperm(sw, so, sg+ss, cells_);
827 
828  Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
829  ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
830  const ADB& po = state.canonical_phase_pressures[ Oil ];
831  const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_)
832  * solvent_props_.pressureMiscibilityFunction(po, cells_);
833 
834  const ADB sn = ss + so + sg;
835 
836  // adjust endpoints
837  const V sgcr = fluid_.scaledCriticalGasSaturations(cells_);
838  const V sogcr = fluid_.scaledCriticalOilinGasSaturations(cells_);
839  const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
840  const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
841 
842  const V ones = V::Constant(nc, 1.0);
843  ADB sor = misc * sorwmis + (ones - misc) * sogcr;
844  ADB sgc = misc * sgcwmis + (ones - misc) * sgcr;
845 
846  ADB ssg = ss + sg - sgc;
847  const ADB sn_eff = sn - sor - sgc;
848 
849  // avoid negative values
850  Selector<double> negSsg_selector(ssg.value(), Selector<double>::LessZero);
851  ssg = negSsg_selector.select(zero, ssg);
852 
853  // avoid negative value and division on zero
854  Selector<double> zeroSn_selector(sn_eff.value(), Selector<double>::LessEqualZero);
855  const ADB F_totalGas = zeroSn_selector.select( zero, ssg / sn_eff);
856 
857  const ADB mkrgt = solvent_props_.miscibleSolventGasRelPermMultiplier(F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
858  const ADB mkro = solvent_props_.miscibleOilRelPermMultiplier(ones - F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_);
859 
860  const V eps = V::Constant(nc, 1e-7);
861  Selector<double> noOil_selector(so.value()-eps, Selector<double>::LessEqualZero);
862  relperm[Gas] = (ones - misc) * relperm[Gas] + misc * mkrgt;
863  relperm[Oil] = noOil_selector.select(relperm[Oil], (ones - misc) * relperm[Oil] + misc * mkro);
864  return relperm;
865  } else {
866  return fluid_.relperm(sw, so, sg+ss, cells_);
867  }
868  } else {
869  return fluid_.relperm(sw, so, sg, cells_);
870  }
871 
872  }
873 
874 
875  template <class Grid>
876  void
877  BlackoilSolventModel<Grid>::computeEffectiveProperties(const SolutionState& state)
878  {
879  // Viscosity
880  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
881  const int np = fluid_.numPhases();
882  const int nc = Opm::UgGridHelpers::numCells(grid_);
883  const ADB zero = ADB::constant(V::Zero(nc));
884 
885  const ADB& pw = state.canonical_phase_pressures[pu.phase_pos[Water]];
886  const ADB& po = state.canonical_phase_pressures[pu.phase_pos[Oil]];
887  const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
888  const std::vector<PhasePresence>& cond = phaseCondition();
889 
890  const ADB mu_w = fluid_.muWat(pw, state.temperature, cells_);
891  const ADB mu_o = fluid_.muOil(po, state.temperature, state.rs, cond, cells_);
892  const ADB mu_g = fluid_.muGas(pg, state.temperature, state.rv, cond, cells_);
893  const ADB mu_s = solvent_props_.muSolvent(pg,cells_);
894  std::vector<ADB> viscosity(np + 1, ADB::null());
895  viscosity[pu.phase_pos[Oil]] = mu_o;
896  viscosity[pu.phase_pos[Gas]] = mu_g;
897  viscosity[pu.phase_pos[Water]] = mu_w;
898  viscosity[solvent_pos_] = mu_s;
899 
900  // Density
901  const ADB bw = fluid_.bWat(pw, state.temperature, cells_);
902  const ADB bo = fluid_.bOil(po, state.temperature, state.rs, cond, cells_);
903  const ADB bg = fluid_.bGas(pg, state.temperature, state.rv, cond, cells_);
904  const ADB bs = solvent_props_.bSolvent(pg, cells_);
905 
906  std::vector<ADB> density(np + 1, ADB::null());
907  density[pu.phase_pos[Oil]] = fluidDensity(Oil, bo, state.rs, state.rv);
908  density[pu.phase_pos[Gas]] = fluidDensity(Gas, bg, state.rs, state.rv);
909  density[pu.phase_pos[Water]] = fluidDensity(Water, bw, state.rs, state.rv);
910  density[solvent_pos_] = fluidDensity(Solvent, bs, state.rs, state.rv);
911 
912  // Effective saturations
913  const ADB& ss = state.solvent_saturation;
914  const ADB& so = state.saturation[ pu.phase_pos[ Oil ] ];
915  const ADB& sg = (active_[ Gas ]
916  ? state.saturation[ pu.phase_pos[ Gas ] ]
917  : zero);
918  const ADB& sw = (active_[ Water ]
919  ? state.saturation[ pu.phase_pos[ Water ] ]
920  : zero);
921 
922  const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_);
923  const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_);
924 
925  std::vector<ADB> effective_saturations (np + 1, ADB::null());
926  effective_saturations[pu.phase_pos[Oil]] = so - sorwmis;
927  effective_saturations[pu.phase_pos[Gas]] = sg - sgcwmis;
928  effective_saturations[pu.phase_pos[Water]] = sw;
929  effective_saturations[solvent_pos_] = ss - sgcwmis;
930 
931  // Compute effective viscosities and densities
932  computeToddLongstaffMixing(viscosity, density, effective_saturations, po, pu);
933 
934  // compute the volume factors from the densities
935  const ADB b_eff_o = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs);
936  const ADB b_eff_g = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv);
937  const ADB b_eff_s = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_);
938 
939  // account for pressure effects and store the computed volume factors and viscosities
940  const V ones = V::Constant(nc, 1.0);
941  const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
942 
943  b_eff_[pu.phase_pos[ Oil ]] = pmisc * b_eff_o + (ones - pmisc) * bo;
944  b_eff_[pu.phase_pos[ Gas ]] = pmisc * b_eff_g + (ones - pmisc) * bg;
945  b_eff_[solvent_pos_] = pmisc * b_eff_s + (ones - pmisc) * bs;
946 
947  // keep the mu*b interpolation
948  mu_eff_[pu.phase_pos[ Oil ]] = b_eff_[pu.phase_pos[ Oil ]] / (pmisc * b_eff_o / viscosity[pu.phase_pos[ Oil ]] + (ones - pmisc) * bo / mu_o);
949  mu_eff_[pu.phase_pos[ Gas ]] = b_eff_[pu.phase_pos[ Gas ]] / (pmisc * b_eff_g / viscosity[pu.phase_pos[ Gas ]] + (ones - pmisc) * bg / mu_g);
950  mu_eff_[solvent_pos_] = b_eff_[solvent_pos_] / (pmisc * b_eff_s / viscosity[solvent_pos_] + (ones - pmisc) * bs / mu_s);
951 
952  // for water the pure values are used
953  mu_eff_[pu.phase_pos[ Water ]] = mu_w;
954  b_eff_[pu.phase_pos[ Water ]] = bw;
955  }
956 
957  template <class Grid>
958  void
959  BlackoilSolventModel<Grid>::computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const ADB po, const Opm::PhaseUsage pu)
960  {
961  const int nc = cells_.size();
962  const V ones = V::Constant(nc, 1.0);
963  const ADB zero = ADB::constant(V::Zero(nc));
964 
965  // Calculation of effective saturations
966  ADB so_eff = saturations[pu.phase_pos[ Oil ]];
967  ADB sg_eff = saturations[pu.phase_pos[ Gas ]];
968  ADB ss_eff = saturations[solvent_pos_];
969 
970  // Avoid negative values
971  Selector<double> negative_selectorSo(so_eff.value(), Selector<double>::LessZero);
972  Selector<double> negative_selectorSg(sg_eff.value(), Selector<double>::LessZero);
973  Selector<double> negative_selectorSs(ss_eff.value(), Selector<double>::LessZero);
974  so_eff = negative_selectorSo.select(zero, so_eff);
975  sg_eff = negative_selectorSg.select(zero, sg_eff);
976  ss_eff = negative_selectorSs.select(zero, ss_eff);
977 
978  // Make copy of the pure viscosities
979  const ADB mu_o = viscosity[pu.phase_pos[ Oil ]];
980  const ADB mu_g = viscosity[pu.phase_pos[ Gas ]];
981  const ADB mu_s = viscosity[solvent_pos_];
982 
983  const ADB sn_eff = so_eff + sg_eff + ss_eff;
984  const ADB sos_eff = so_eff + ss_eff;
985  const ADB ssg_eff = ss_eff + sg_eff;
986 
987  // Avoid division by zero
988  Selector<double> zero_selectorSos(sos_eff.value(), Selector<double>::Zero);
989  Selector<double> zero_selectorSsg(ssg_eff.value(), Selector<double>::Zero);
990  Selector<double> zero_selectorSn(sn_eff.value(), Selector<double>::Zero);
991 
992  const ADB mu_s_pow = pow(mu_s, 0.25);
993  const ADB mu_o_pow = pow(mu_o, 0.25);
994  const ADB mu_g_pow = pow(mu_g, 0.25);
995 
996  const ADB mu_mos = zero_selectorSos.select(mu_o, mu_o * mu_s / pow( ( (so_eff / sos_eff) * mu_s_pow) + ( (ss_eff / sos_eff) * mu_o_pow) , 4.0));
997  const ADB mu_msg = zero_selectorSsg.select(mu_g, mu_g * mu_s / pow( ( (sg_eff / ssg_eff) * mu_s_pow) + ( (ss_eff / ssg_eff) * mu_g_pow) , 4.0));
998  const ADB mu_m = zero_selectorSn.select(mu_s, mu_o * mu_s * mu_g / pow( ( (so_eff / sn_eff) * mu_s_pow * mu_g_pow)
999  + ( (ss_eff / sn_eff) * mu_o_pow * mu_g_pow) + ( (sg_eff / sn_eff) * mu_s_pow * mu_o_pow), 4.0));
1000  // Mixing parameter for viscosity
1001  // The pressureMixingParameter represent the miscibility of the solvent while the mixingParameterViscosity the effect of the porous media.
1002  // The pressureMixingParameter is not implemented in ecl100.
1003  const ADB mix_param_mu = solvent_props_.mixingParameterViscosity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
1004 
1005  Selector<double> zero_selectorSs(ss_eff.value(), Selector<double>::Zero);
1006  // Update viscosities, use pure values if solvent saturation is zero
1007  viscosity[pu.phase_pos[ Oil ]] = zero_selectorSs.select(mu_o, pow(mu_o,ones - mix_param_mu) * pow(mu_mos, mix_param_mu));
1008  viscosity[pu.phase_pos[ Gas ]] = zero_selectorSs.select(mu_g, pow(mu_g,ones - mix_param_mu) * pow(mu_msg, mix_param_mu));
1009  viscosity[solvent_pos_] = zero_selectorSs.select(mu_s, pow(mu_s,ones - mix_param_mu) * pow(mu_m, mix_param_mu));
1010 
1011  // Density
1012  ADB& rho_o = density[pu.phase_pos[ Oil ]];
1013  ADB& rho_g = density[pu.phase_pos[ Gas ]];
1014  ADB& rho_s = density[solvent_pos_];
1015 
1016  // mixing parameter for density
1017  const ADB mix_param_rho = solvent_props_.mixingParameterDensity(cells_) * solvent_props_.pressureMixingParameter(po, cells_);
1018 
1019  // compute effective viscosities for density calculations. These have to
1020  // be recomputed as a different mixing parameter may be used.
1021  const ADB mu_o_eff = pow(mu_o,ones - mix_param_rho) * pow(mu_mos, mix_param_rho);
1022  const ADB mu_g_eff = pow(mu_g,ones - mix_param_rho) * pow(mu_msg, mix_param_rho);
1023  const ADB mu_s_eff = pow(mu_s,ones - mix_param_rho) * pow(mu_m, mix_param_rho);
1024 
1025  const ADB sog_eff = so_eff + sg_eff;
1026  // Avoid division by zero
1027  Selector<double> zero_selectorSog_eff(sog_eff.value(), Selector<double>::Zero);
1028  const ADB sof = zero_selectorSog_eff.select(zero , so_eff / sog_eff);
1029  const ADB sgf = zero_selectorSog_eff.select(zero , sg_eff / sog_eff);
1030 
1031  // Effective densities
1032  const ADB mu_sog_pow = mu_s_pow * ( (sgf * mu_o_pow) + (sof * mu_g_pow) );
1033  const ADB mu_o_eff_pow = pow(mu_o_eff, 0.25);
1034  const ADB mu_g_eff_pow = pow(mu_g_eff, 0.25);
1035  const ADB mu_s_eff_pow = pow(mu_s_eff, 0.25);
1036  const ADB sfraction_oe = (mu_o_pow * (mu_o_eff_pow - mu_s_pow)) / (mu_o_eff_pow * (mu_o_pow - mu_s_pow));
1037  const ADB sfraction_ge = (mu_g_pow * (mu_s_pow - mu_g_eff_pow)) / (mu_g_eff_pow * (mu_s_pow - mu_g_pow));
1038  const ADB sfraction_se = (mu_sog_pow - ( mu_o_pow * mu_g_pow * mu_s_pow / mu_s_eff_pow) ) / ( mu_sog_pow - (mu_o_pow * mu_g_pow));
1039  const ADB rho_o_eff = (rho_o * sfraction_oe) + (rho_s * (ones - sfraction_oe));
1040  const ADB rho_g_eff = (rho_g * sfraction_ge) + (rho_s * (ones - sfraction_ge));
1041  const ADB rho_s_eff = (rho_s * sfraction_se) + (rho_g * sgf * (ones - sfraction_se)) + (rho_o * sof * (ones - sfraction_se));
1042 
1043  // Avoid division by zero for equal mobilities. For equal mobilities the effecitive density is calculated
1044  // based on the saturation fraction directly.
1045  Selector<double> unitGasSolventMobilityRatio_selector(mu_s.value() - mu_g.value(), Selector<double>::Zero);
1046  Selector<double> unitOilSolventMobilityRatio_selector(mu_s.value() - mu_o.value(), Selector<double>::Zero);
1047 
1048  // Effective densities when the mobilities are equal
1049  const ADB rho_m = zero_selectorSn.select(zero, (rho_o * so_eff / sn_eff) + (rho_g * sg_eff / sn_eff) + (rho_s * ss_eff / sn_eff));
1050  const ADB rho_o_eff_simple = ((ones - mix_param_rho) * rho_o) + (mix_param_rho * rho_m);
1051  const ADB rho_g_eff_simple = ((ones - mix_param_rho) * rho_g) + (mix_param_rho * rho_m);
1052  //const ADB rho_s_eff_simple = ((ones - mix_param_rho) * rho_s) + (mix_param_rho * rho_m);
1053 
1054  // Update densities, use pure values if solvent saturation is zero
1055  rho_o = zero_selectorSs.select(rho_o, unitOilSolventMobilityRatio_selector.select(rho_o_eff_simple, rho_o_eff) );
1056  rho_g = zero_selectorSs.select(rho_g, unitGasSolventMobilityRatio_selector.select(rho_g_eff_simple, rho_g_eff) );
1057  rho_s = zero_selectorSs.select(rho_s, rho_s_eff);
1058 
1059 
1060  }
1061 
1062 
1063  template <class Grid>
1064  std::vector<ADB>
1065  BlackoilSolventModel<Grid>::
1066  computePressures(const ADB& po,
1067  const ADB& sw,
1068  const ADB& so,
1069  const ADB& sg,
1070  const ADB& ss) const
1071  {
1072  std::vector<ADB> pressures = Base::computePressures(po, sw, so, sg);
1073 
1074  if (has_solvent_) {
1075  // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss)
1076  std::vector<ADB> pressures_imisc = Base::computePressures(po, sw, so, sg + ss);
1077 
1078  // Pressure effects on capillary pressure miscibility
1079  const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
1080  // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent
1081  // to changing the gas pressure directly.
1082  const int nc = cells_.size();
1083  const V ones = V::Constant(nc, 1.0);
1084  pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas]));
1085  }
1086  return pressures;
1087  }
1088 
1089 
1090 
1091 
1092  template <class Grid>
1093  std::vector<std::vector<double> >
1094  BlackoilSolventModel<Grid>::
1095  computeFluidInPlace(const ReservoirState& x,
1096  const std::vector<int>& fipnum)
1097  {
1098  if (has_solvent_ && is_miscible_ && b_eff_[0].size() == 0) {
1099  // A hack to avoid trouble for initial fluid in place, due to usage of b_eff_.
1100  WellState xw, xwdummy;
1101  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
1102  xw.init(&wells(), x, xwdummy, pu);
1103  SolutionState solstate = variableState(x, xw);
1104  computeEffectiveProperties(solstate);
1105  }
1106  return Base::computeFluidInPlace(x, fipnum);
1107  }
1108 
1109 
1110 
1111 }
1112 
1113 
1114 #endif // OPM_BLACKOILSOLVENT_IMPL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
BlackoilSolventModel(const typename Base::ModelParameters &param, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const SolventPropsAdFromDeck &solvent_props, const StandardWellsSolvent &well_model, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output, const bool has_solvent, const bool is_miscible)
Construct the model.
Definition: BlackoilSolventModel_impl.hpp:73
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
Class for handling the standard well model for solvent model.
Definition: StandardWellsSolvent.hpp:32
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilSolventModel_impl.hpp:389
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:364
AutoDiffBlock< Scalar > pow(const AutoDiffBlock< Scalar > &base, const double exponent)
Computes the value of base raised to the power of exponent.
Definition: AutoDiffBlock.hpp:636
std::vector< ADB > material_balance_eq
The material_balance_eq vector has one element for each active phase, each of which has size equal to...
Definition: LinearisedBlackoilResidual.hpp:54
Definition: GridHelpers.cpp:27
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
WellModel & wellModel()
return the WellModel object
Definition: BlackoilModelBase.hpp:268
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
Definition: AutoDiffHelpers.hpp:623
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
Definition: SolventPropsAdFromDeck.hpp:37
ADB select(const ADB &x1, const ADB &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:415
Eigen::Array< double, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
int numPhases() const
Definition: BlackoilPropsAdFromDeck.cpp:225
A model implementation for three-phase black oil with one extra component.
Definition: BlackoilSolventModel.hpp:38
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
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