All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
BlackoilPolymerModel_impl.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2014, 2015 Statoil ASA.
5  Copyright 2015 NTNU
6  Copyright 2015 IRIS AS
7 
8  This file is part of the Open Porous Media project (OPM).
9 
10  OPM is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  OPM is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with OPM. If not, see <http://www.gnu.org/licenses/>.
22 */
23 
24 #ifndef OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
25 #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
26 
27 #include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
28 
29 #include <opm/autodiff/AutoDiffBlock.hpp>
30 #include <opm/autodiff/AutoDiffHelpers.hpp>
31 #include <opm/autodiff/GridHelpers.hpp>
32 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
33 #include <opm/autodiff/GeoProps.hpp>
34 #include <opm/autodiff/WellDensitySegmented.hpp>
35 
36 #include <opm/core/grid.h>
37 #include <opm/core/linalg/LinearSolverInterface.hpp>
38 #include <opm/core/linalg/ParallelIstlInformation.hpp>
39 #include <opm/core/props/rock/RockCompressibility.hpp>
40 #include <opm/common/ErrorMacros.hpp>
41 #include <opm/common/Exceptions.hpp>
42 #include <opm/parser/eclipse/Units/Units.hpp>
43 #include <opm/core/well_controls.h>
44 #include <opm/core/utility/parameters/ParameterGroup.hpp>
45 
46 #include <cassert>
47 #include <cmath>
48 #include <iostream>
49 #include <iomanip>
50 #include <limits>
51 
52 namespace Opm {
53 
54 
55 
56  namespace detail {
57 
58  template <class PU>
59  int polymerPos(const PU& pu)
60  {
61  const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
62  int pos = 0;
63  for (int phase = 0; phase < maxnp; ++phase) {
64  if (pu.phase_used[phase]) {
65  pos++;
66  }
67  }
68 
69  return pos;
70  }
71 
72  } // namespace detail
73 
74 
75 
76  template <class Grid>
77  BlackoilPolymerModel<Grid>::BlackoilPolymerModel(const typename Base::ModelParameters& param,
78  const Grid& grid,
79  const BlackoilPropsAdFromDeck& fluid,
80  const DerivedGeology& geo,
81  const RockCompressibility* rock_comp_props,
82  const PolymerPropsAd& polymer_props_ad,
83  const StandardWells& well_model,
84  const NewtonIterationBlackoilInterface& linsolver,
85  std::shared_ptr< const EclipseState > eclipse_state,
86  const bool has_disgas,
87  const bool has_vapoil,
88  const bool has_polymer,
89  const bool has_plyshlog,
90  const bool has_shrate,
91  const std::vector<double>& wells_rep_radius,
92  const std::vector<double>& wells_perf_length,
93  const std::vector<double>& wells_bore_diameter,
94  const bool terminal_output)
95  : Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver, eclipse_state,
96  has_disgas, has_vapoil, terminal_output),
97  polymer_props_ad_(polymer_props_ad),
98  has_polymer_(has_polymer),
99  has_plyshlog_(has_plyshlog),
100  has_shrate_(has_shrate),
101  poly_pos_(detail::polymerPos(fluid.phaseUsage())),
102  wells_rep_radius_(wells_rep_radius),
103  wells_perf_length_(wells_perf_length),
104  wells_bore_diameter_(wells_bore_diameter)
105  {
106  if (has_polymer_) {
107  if (!active_[Water]) {
108  OPM_THROW(std::logic_error, "Polymer must solved in water!\n");
109  }
110  residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169); // use the same as the water phase
111  // If deck has polymer, residual_ should contain polymer equation.
112  sd_.rq.resize(fluid_.numPhases() + 1);
113  residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
114  Base::material_name_.push_back("Polymer");
115  assert(poly_pos_ == fluid_.numPhases());
116  }
117  }
118 
119 
120 
121  template <class Grid>
122  void
125  const ReservoirState& reservoir_state,
126  const WellState& well_state)
127  {
128  Base::prepareStep(timer, reservoir_state, well_state);
129  auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
130  // Initial max concentration of this time step from PolymerBlackoilState.
131 
132  cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
133  }
134 
135 
136 
137 
138  template <class Grid>
139  void
142  ReservoirState& reservoir_state,
143  WellState& /* well_state */)
144  {
145  computeCmax(reservoir_state);
146  }
147 
148 
149 
150 
151 
152  template <class Grid>
153  void
154  BlackoilPolymerModel<Grid>::makeConstantState(SolutionState& state) const
155  {
156  Base::makeConstantState(state);
157  state.concentration = ADB::constant(state.concentration.value());
158  }
159 
160 
161 
162 
163 
164  template <class Grid>
165  std::vector<V>
166  BlackoilPolymerModel<Grid>::variableStateInitials(const ReservoirState& x,
167  const WellState& xw) const
168  {
169  std::vector<V> vars0 = Base::variableStateInitials(x, xw);
170  assert(int(vars0.size()) == fluid_.numPhases() + 2);
171 
172  // Initial polymer concentration.
173  if (has_polymer_) {
174  const auto& concentration = x.getCellData( x.CONCENTRATION );
175  assert (not concentration.empty());
176  const int nc = concentration.size();
177  const V c = Eigen::Map<const V>(concentration.data() , nc);
178  // Concentration belongs after other reservoir vars but before well vars.
179  auto concentration_pos = vars0.begin() + fluid_.numPhases();
180  assert(concentration_pos == vars0.end() - 2);
181  vars0.insert(concentration_pos, c);
182  }
183  return vars0;
184  }
185 
186 
187 
188 
189 
190  template <class Grid>
191  std::vector<int>
192  BlackoilPolymerModel<Grid>::variableStateIndices() const
193  {
194  std::vector<int> ind = Base::variableStateIndices();
195  assert(ind.size() == 5);
196  if (has_polymer_) {
197  ind.resize(6);
198  // Concentration belongs after other reservoir vars but before well vars.
199  ind[Concentration] = fluid_.numPhases();
200  // Concentration is pushing back the well vars.
201  ++ind[Qs];
202  ++ind[Bhp];
203  }
204  return ind;
205  }
206 
207 
208 
209 
210  template <class Grid>
211  typename BlackoilPolymerModel<Grid>::SolutionState
212  BlackoilPolymerModel<Grid>::variableStateExtractVars(const ReservoirState& x,
213  const std::vector<int>& indices,
214  std::vector<ADB>& vars) const
215  {
216  SolutionState state = Base::variableStateExtractVars(x, indices, vars);
217  if (has_polymer_) {
218  state.concentration = std::move(vars[indices[Concentration]]);
219  }
220  return state;
221  }
222 
223 
224 
225 
226 
227  template <class Grid>
228  void
229  BlackoilPolymerModel<Grid>::computeAccum(const SolutionState& state,
230  const int aix )
231  {
232  Base::computeAccum(state, aix);
233 
234  // Compute accumulation of polymer equation only if needed.
235  if (has_polymer_) {
236  const ADB& press = state.pressure;
237  const std::vector<ADB>& sat = state.saturation;
238  const ADB& c = state.concentration;
239  const ADB pv_mult = poroMult(press); // also computed in Base::computeAccum, could be optimized.
240  const Opm::PhaseUsage& pu = fluid_.phaseUsage();
241  // compute polymer properties.
242  const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
243  const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
244  const double rho_rock = polymer_props_ad_.rockDensity();
245  const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], AutoDiffGrid::numCells(grid_));
246  const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
247  // Compute polymer accumulation term.
248  sd_.rq[poly_pos_].accum[aix] = pv_mult * sd_.rq[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
249  + pv_mult * rho_rock * (1. - phi) / phi * ads;
250  }
251 
252  }
253 
254 
255 
256 
257 
258 
259  template <class Grid>
260  void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
261  {
262  auto& max_concentration = state.getCellData( state.CMAX );
263  const auto& concentration = state.getCellData( state.CONCENTRATION );
264  std::transform( max_concentration.begin() ,
265  max_concentration.end() ,
266  concentration.begin() ,
267  max_concentration.begin() ,
268  [](double c_max , double c) { return std::max( c_max , c ); });
269  }
270 
271 
272 
273 
274 
275  template <class Grid>
276  void
277  BlackoilPolymerModel<Grid>::
278  assembleMassBalanceEq(const SolutionState& state)
279  {
280  // Base::assembleMassBalanceEq(state);
281 
282  // Compute b_p and the accumulation term b_p*s_p for each phase,
283  // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
284  // These quantities are stored in sd_.rq[phase].accum[1].
285  // The corresponding accumulation terms from the start of
286  // the timestep (b^0_p*s^0_p etc.) were already computed
287  // on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
288  computeAccum(state, 1);
289 
290  // Set up the common parts of the mass balance equations
291  // for each active phase.
292  const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
293  {
294  const std::vector<ADB> kr = computeRelPerm(state);
295  for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
296  sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
297  }
298  }
299 
300 
301  if (has_plyshlog_) {
302  std::vector<double> water_vel;
303  std::vector<double> visc_mult;
304 
305  computeWaterShearVelocityFaces(transi, state.canonical_phase_pressures, state, water_vel, visc_mult);
306  if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
307  // std::cerr << " failed in calculating the shear-multiplier " << std::endl;
308  OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
309  }
310  }
311 
312  for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
313  const std::vector<PhasePresence>& cond = phaseCondition();
314  sd_.rq[phaseIdx].mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
315  sd_.rq[phaseIdx].rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
316  computeMassFlux(phaseIdx, transi, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
317 
318  residual_.material_balance_eq[ phaseIdx ] =
319  pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
320  + ops_.div*sd_.rq[phaseIdx].mflux;
321  }
322 
323  // -------- Extra (optional) rs and rv contributions to the mass balance equations --------
324 
325  // Add the extra (flux) terms to the mass balance equations
326  // From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
327  // The extra terms in the accumulation part of the equation are already handled.
328  if (active_[ Oil ] && active_[ Gas ]) {
329  const int po = fluid_.phaseUsage().phase_pos[ Oil ];
330  const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
331 
332  const UpwindSelector<double> upwindOil(grid_, ops_,
333  sd_.rq[po].dh.value());
334  const ADB rs_face = upwindOil.select(state.rs);
335 
336  const UpwindSelector<double> upwindGas(grid_, ops_,
337  sd_.rq[pg].dh.value());
338  const ADB rv_face = upwindGas.select(state.rv);
339 
340  residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
341  residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
342 
343  // OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
344 
345  }
346 
347  // Add polymer equation.
348  if (has_polymer_) {
349  residual_.material_balance_eq[poly_pos_] = pvdt_ * (sd_.rq[poly_pos_].accum[1] - sd_.rq[poly_pos_].accum[0])
350  + ops_.div*sd_.rq[poly_pos_].mflux;
351  }
352 
353 
354  if (param_.update_equations_scaling_) {
355  updateEquationsScaling();
356  }
357 
358  }
359 
360 
361 
362 
363 
364  template <class Grid>
365  void BlackoilPolymerModel<Grid>::updateEquationsScaling()
366  {
367  Base::updateEquationsScaling();
368  if (has_polymer_) {
369  const int water_pos = fluid_.phaseUsage().phase_pos[Water];
370  residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos];
371  }
372  }
373 
374 
375 
376 
377 
378  template <class Grid>
379  void BlackoilPolymerModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
380  const SolutionState& state,
381  WellState& xw)
382 
383  {
384  Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
385 
386  // Add well contributions to polymer mass balance equation
387  if (has_polymer_) {
388  const ADB mc = computeMc(state);
389  const int nc = xw.polymerInflow().size();
390  const V polyin = Eigen::Map<const V>(xw.polymerInflow().data(), nc);
391  const int nperf = wells().well_connpos[wells().number_of_wells];
392  const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
393  const V poly_in_perf = subset(polyin, well_cells);
394  const V poly_mc_perf = subset(mc.value(), well_cells);
395  const ADB& cq_s_water = cq_s[fluid_.phaseUsage().phase_pos[Water]];
396  Selector<double> injector_selector(cq_s_water.value());
397  const V poly_perf = injector_selector.select(poly_in_perf, poly_mc_perf);
398  const ADB cq_s_poly = cq_s_water * poly_perf;
399  residual_.material_balance_eq[poly_pos_] -= superset(cq_s_poly, well_cells, nc);
400  }
401  }
402 
403 
404 
405 
406 
407 
408  template <class Grid>
410  ReservoirState& reservoir_state,
411  WellState& well_state)
412  {
413  if (has_polymer_) {
414  // Extract concentration change.
415  const int np = fluid_.numPhases();
416  const int nc = Opm::AutoDiffGrid::numCells(grid_);
417  const V zero = V::Zero(nc);
418  const int concentration_start = nc * np;
419  const V dc = subset(dx, Span(nc, 1, concentration_start));
420 
421  // Create new dx with the dc part deleted.
422  V modified_dx = V::Zero(dx.size() - nc);
423  modified_dx.head(concentration_start) = dx.head(concentration_start);
424  const int tail_len = dx.size() - concentration_start - nc;
425  modified_dx.tail(tail_len) = dx.tail(tail_len);
426 
427  // Call base version.
428  Base::updateState(modified_dx, reservoir_state, well_state);
429 
430  {
431  auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
432  // Update concentration.
433  const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
434  const V c = (c_old - dc).max(zero);
435  std::copy(&c[0], &c[0] + nc, concentration.begin());
436  }
437  } else {
438  // Just forward call to base version.
439  Base::updateState(dx, reservoir_state, well_state);
440  }
441  }
442 
443 
444 
445 
446 
447  template <class Grid>
448  void
450  const V& transi,
451  const ADB& kr ,
452  const ADB& mu ,
453  const ADB& rho ,
454  const ADB& phasePressure,
455  const SolutionState& state)
456  {
457  Base::computeMassFlux(actph, transi, kr, mu, rho, phasePressure, state);
458 
459  // Polymer treatment.
460  const int canonicalPhaseIdx = canph_[ actph ];
461  if (canonicalPhaseIdx == Water) {
462  if (has_polymer_) {
463  const ADB tr_mult = transMult(state.pressure);
464  const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
465  const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
466  const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
467  // Reduce mobility of water phase by relperm reduction and effective viscosity increase.
468  sd_.rq[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
469  // Compute polymer mobility.
470  const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value());
471  sd_.rq[poly_pos_].mob = tr_mult * krw_eff * state.concentration * inv_poly_eff_visc;
472  sd_.rq[poly_pos_].b = sd_.rq[actph].b;
473  sd_.rq[poly_pos_].dh = sd_.rq[actph].dh;
474  UpwindSelector<double> upwind(grid_, ops_, sd_.rq[poly_pos_].dh.value());
475  // Compute polymer flux.
476  sd_.rq[poly_pos_].mflux = upwind.select(sd_.rq[poly_pos_].b * sd_.rq[poly_pos_].mob) * (transi * sd_.rq[poly_pos_].dh);
477  // Must recompute water flux since we have to use modified mobilities.
478  sd_.rq[ actph ].mflux = upwind.select(sd_.rq[actph].b * sd_.rq[actph].mob) * (transi * sd_.rq[actph].dh);
479 
480  // applying the shear-thinning factors
481  if (has_plyshlog_) {
482  V shear_mult_faces_v = Eigen::Map<V>(shear_mult_faces_.data(), shear_mult_faces_.size());
483  ADB shear_mult_faces_adb = ADB::constant(shear_mult_faces_v);
484  sd_.rq[poly_pos_].mflux = sd_.rq[poly_pos_].mflux / shear_mult_faces_adb;
485  sd_.rq[actph].mflux = sd_.rq[actph].mflux / shear_mult_faces_adb;
486  }
487  }
488  }
489  }
490 
491 
492 
493 
494  template <class Grid>
495  SimulatorReport
496  BlackoilPolymerModel<Grid>::assemble(const ReservoirState& reservoir_state,
497  WellState& well_state,
498  const bool initial_assembly)
499  {
500  using namespace Opm::AutoDiffGrid;
501 
502  SimulatorReport report;
503 
504  // Possibly switch well controls and updating well state to
505  // get reasonable initial conditions for the wells
506  // updateWellControls(well_state);
507  wellModel().updateWellControls(well_state);
508 
509  // Create the primary variables.
510  SolutionState state = variableState(reservoir_state, well_state);
511 
512  if (initial_assembly) {
513  // Create the (constant, derivativeless) initial state.
514  SolutionState state0 = state;
515  makeConstantState(state0);
516  // Compute initial accumulation contributions
517  // and well connection pressures.
518  computeAccum(state0, 0);
519  // computeWellConnectionPressures(state0, well_state);
520  wellModel().computeWellConnectionPressures(state0, well_state);
521  }
522 
523  // OPM_AD_DISKVAL(state.pressure);
524  // OPM_AD_DISKVAL(state.saturation[0]);
525  // OPM_AD_DISKVAL(state.saturation[1]);
526  // OPM_AD_DISKVAL(state.saturation[2]);
527  // OPM_AD_DISKVAL(state.rs);
528  // OPM_AD_DISKVAL(state.rv);
529  // OPM_AD_DISKVAL(state.qs);
530  // OPM_AD_DISKVAL(state.bhp);
531 
532  // -------- Mass balance equations --------
533  assembleMassBalanceEq(state);
534  // -------- Well equations ----------
535  if ( ! wellsActive() ) {
536  return report;
537  }
538 
539  std::vector<ADB> mob_perfcells;
540  std::vector<ADB> b_perfcells;
541  wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
542 
543  // updating the the injection mobility related to polymer injection when necessary
544  // only the mobility of water phase is updated
545  computeInjectionMobility(state, mob_perfcells);
546 
547  if (param_.solve_welleq_initially_ && initial_assembly) {
548  // solve the well equations as a pre-processing step
549  Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
550  }
551 
552  V aliveWells;
553  std::vector<ADB> cq_s;
554 
555  wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
556 
557  if (has_plyshlog_) {
558  std::vector<double> water_vel_wells;
559  std::vector<double> visc_mult_wells;
560 
561  const int water_pos = fluid_.phaseUsage().phase_pos[Water];
562  computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells);
563 
564  if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) {
565  OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells ");
566  }
567 
568  // applying the shear-thinning to the water phase
569  V shear_mult_wells_v = Eigen::Map<V>(shear_mult_wells_.data(), shear_mult_wells_.size());
570  ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v);
571  mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
572  wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
573  }
574 
575  wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
576  wellModel().addWellFluxEq(cq_s, state, residual_);
577  addWellContributionToMassBalanceEq(cq_s, state, well_state);
578  wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
579 
580  report.converged = true;
581  return report;
582  }
583 
584 
585 
586 
587  template <class Grid>
588  ADB
589  BlackoilPolymerModel<Grid>::computeMc(const SolutionState& state) const
590  {
591  return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
592  }
593 
594 
595 
596 
597  template<class Grid>
598  void
600  const std::vector<ADB>& phasePressure, const SolutionState& state,
601  std::vector<double>& water_vel, std::vector<double>& visc_mult)
602  {
603 
604  const int phase = fluid_.phaseUsage().phase_pos[Water]; // water position
605 
606  const int canonicalPhaseIdx = canph_[phase];
607 
608  const std::vector<PhasePresence> cond = phaseCondition();
609 
610  const ADB tr_mult = transMult(state.pressure);
611  const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
612  sd_.rq[phase].mob = tr_mult * sd_.rq[phase].kr / mu;
613 
614  // compute gravity potensial using the face average as in eclipse and MRST
615  const ADB rho = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
616  const ADB rhoavg = ops_.caver * rho;
617  sd_.rq[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
618  if (use_threshold_pressure_) {
619  applyThresholdPressures(sd_.rq[ phase ].dh);
620  }
621 
622  const ADB& b = sd_.rq[ phase ].b;
623  const ADB& mob = sd_.rq[ phase ].mob;
624  const ADB& dh = sd_.rq[ phase ].dh;
625  UpwindSelector<double> upwind(grid_, ops_, dh.value());
626 
627  const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
628  ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
629  cmax,
630  sd_.rq[phase].kr);
631  ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
632  sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
633 
634  const V& polymer_conc = state.concentration.value();
635  V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
636  V visc_mult_faces = upwind.select(visc_mult_cells);
637 
638  size_t nface = visc_mult_faces.size();
639  visc_mult.resize(nface);
640  std::copy(visc_mult_faces.data(), visc_mult_faces.data() + nface, visc_mult.begin());
641 
642  sd_.rq[ phase ].mflux = (transi * upwind.select(b * mob)) * dh;
643 
644 
645  const auto& b_faces_adb = upwind.select(b);
646  std::vector<double> b_faces(b_faces_adb.value().data(), b_faces_adb.value().data() + b_faces_adb.size());
647 
648  const auto& internal_faces = ops_.internal_faces;
649 
650  std::vector<double> internal_face_areas;
651  internal_face_areas.resize(internal_faces.size());
652 
653  for (int i = 0; i < internal_faces.size(); ++i) {
654  internal_face_areas[i] = grid_.face_areas[internal_faces[i]];
655  }
656 
657  const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
658  const ADB phiavg_adb = ops_.caver * phi;
659 
660  std::vector<double> phiavg(phiavg_adb.value().data(), phiavg_adb.value().data() + phiavg_adb.size());
661 
662  water_vel.resize(nface);
663  std::copy(sd_.rq[0].mflux.value().data(), sd_.rq[0].mflux.value().data() + nface, water_vel.begin());
664 
665  for (size_t i = 0; i < nface; ++i) {
666  water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]);
667  }
668 
669  // for SHRATE keyword treatment
670  if (has_shrate_) {
671 
672  // get the upwind water saturation
673  const Opm::PhaseUsage pu = fluid_.phaseUsage();
674  const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
675  const ADB& sw_upwind_adb = upwind.select(sw);
676  std::vector<double> sw_upwind(sw_upwind_adb.value().data(), sw_upwind_adb.value().data() + sw_upwind_adb.size());
677 
678  // get the absolute permeability for the faces
679  std::vector<double> perm;
680  perm.resize(transi.size());
681 
682  for (int i = 0; i < transi.size(); ++i) {
683  perm[i] = transi[i] / internal_faces[i];
684  }
685 
686  // get the upwind krw_eff
687  const ADB& krw_adb = upwind.select(krw_eff);
688  std::vector<double> krw_upwind(krw_adb.value().data(), krw_adb.value().data() + krw_adb.size());
689 
690  const double& shrate_const = polymer_props_ad_.shrate();
691 
692  const double epsilon = std::numeric_limits<double>::epsilon();
693  // std::cout << "espilon is " << epsilon << std::endl;
694  // std::cin.ignore();
695 
696  for (size_t i = 0; i < water_vel.size(); ++i) {
697  // assuming only when upwinding water saturation is not zero
698  // there will be non-zero water velocity
699  if (std::abs(water_vel[i]) < epsilon) {
700  continue;
701  }
702 
703  water_vel[i] *= shrate_const * std::sqrt(phiavg[i] / (perm[i] * sw_upwind[i] * krw_upwind[i]));
704 
705  }
706  }
707 
708  }
709 
710 
711 
712 
713  template<class Grid>
714  void
715  BlackoilPolymerModel<Grid>::computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw,
716  std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells)
717  {
718  if( ! wellsActive() ) return ;
719 
720  const int nw = wells().number_of_wells;
721  const int nperf = wells().well_connpos[nw];
722  const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
723 
724  water_vel_wells.resize(cq_sw.size());
725  std::copy(cq_sw.value().data(), cq_sw.value().data() + cq_sw.size(), water_vel_wells.begin());
726 
727  const V& polymer_conc = state.concentration.value();
728 
729  V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
730  V visc_mult_wells_v = subset(visc_mult_cells, well_cells);
731 
732  visc_mult_wells.resize(visc_mult_wells_v.size());
733  std::copy(visc_mult_wells_v.data(), visc_mult_wells_v.data() + visc_mult_wells_v.size(), visc_mult_wells.begin());
734 
735  const int water_pos = fluid_.phaseUsage().phase_pos[Water];
736  ADB b_perfcells = subset(sd_.rq[water_pos].b, well_cells);
737 
738  const ADB& p_perfcells = subset(state.pressure, well_cells);
739  const V& cdp = wellModel().wellPerforationPressureDiffs();
740  const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
741  // Pressure drawdown (also used to determine direction of flow)
742  const ADB drawdown = p_perfcells - perfpressure;
743 
744  // selects injection perforations
745  V selectInjectingPerforations = V::Zero(nperf);
746  for (int c = 0; c < nperf; ++c) {
747  if (drawdown.value()[c] < 0) {
748  selectInjectingPerforations[c] = 1;
749  }
750  }
751 
752  // for the injection wells
753  for (size_t i = 0; i < well_cells.size(); ++i) {
754  if (xw.polymerInflow()[well_cells[i]] == 0. && selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
755  visc_mult_wells[i] = 1.;
756  }
757  if (selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
758  if (xw.polymerInflow()[well_cells[i]] == 0.) {
759  visc_mult_wells[i] = 1.;
760  } else {
761  // TODO: not tested for this assumption yet
762  const double c_perf = state.concentration.value()[well_cells[i]];
763  visc_mult_wells[i] = polymer_props_ad_.viscMult(c_perf);
764  }
765  }
766  }
767 
768  const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
769  const ADB phi_wells_adb = subset(phi, well_cells);
770 
771  std::vector<double> phi_wells(phi_wells_adb.value().data(), phi_wells_adb.value().data() + phi_wells_adb.size());
772 
773  std::vector<double> b_wells(b_perfcells.value().data(), b_perfcells.value().data() + b_perfcells.size());
774 
775  for (size_t i = 0; i < water_vel_wells.size(); ++i) {
776  water_vel_wells[i] = b_wells[i] * water_vel_wells[i] / (phi_wells[i] * 2. * M_PI * wells_rep_radius_[i] * wells_perf_length_[i]);
777  // TODO: CHECK to make sure this formulation is corectly used. Why muliplied by bW.
778  // Although this formulation works perfectly with the tests compared with other formulations
779  }
780 
781  // for SHRATE treatment
782  if (has_shrate_) {
783  const double& shrate_const = polymer_props_ad_.shrate();
784  for (size_t i = 0; i < water_vel_wells.size(); ++i) {
785  water_vel_wells[i] = shrate_const * water_vel_wells[i] / wells_bore_diameter_[i];
786  }
787  }
788 
789  return;
790  }
791 
792 
793 
794 
795 
796  template<class Grid>
797  void
799  computeInjectionMobility(const SolutionState& state,
800  std::vector<ADB>& mob_perfcells)
801  {
802  // Calculating the mobility for the polymer injection peforations
803  if (has_polymer_ && wellModel().localWellsActive()) {
804  const std::vector<int> well_cells = wellModel().wellOps().well_cells;
805  const int nperf = well_cells.size();
806 
807  // Calculating the drawdown to decide the injection perforation
808  const ADB& p_perfcells = subset(state.pressure, well_cells);
809  const V& cdp = wellModel().wellPerforationPressureDiffs();
810  const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
811  // Pressure drawdown (also used to determine direction of flow)
812  const ADB drawdown = p_perfcells - perfpressure;
813 
814  // Polymer concentration in the perforations
815  const ADB c_perfcells = subset(state.concentration, well_cells);
816 
817  // Distinguishing the injection perforation from other perforation
818  // The value is the location in the well_cell array, not the global index
819  std::vector<int> polymer_inj_cells;
820  std::vector<int> other_well_cells;
821 
822  polymer_inj_cells.reserve(nperf);
823  other_well_cells.reserve(nperf);
824 
825  for (int c = 0; c < nperf; ++c) {
826  // TODO: more tests need to be done for this criterion
827  if (drawdown.value()[c] < 0.0 && c_perfcells.value()[c] > 0.0) {
828  polymer_inj_cells.push_back(c);
829  } else {
830  other_well_cells.push_back(c);
831  }
832  }
833 
834  // there is some polymer injection process going
835  if ( !polymer_inj_cells.empty() ) {
836  // the mobility need to be recalculated for the polymer injection cells
837  const int water_pos = fluid_.phaseUsage().phase_pos[Water];
838  const ADB mu_perfcells = subset(sd_.rq[water_pos].mu, well_cells);
839 
840  const ADB c_poly_inj_cells = subset(c_perfcells, polymer_inj_cells);
841  const ADB mu_poly_inj_cells = subset(mu_perfcells, polymer_inj_cells); // water viscosity
842 
843  const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(c_poly_inj_cells, mu_poly_inj_cells.value());
844  const ADB fully_mixing_visc = polymer_props_ad_.viscMult(c_poly_inj_cells) * mu_poly_inj_cells;
845 
846  // the original mobility for the polymer injection well cells
847  ADB mob_polymer_inj = subset(mob_perfcells[water_pos], polymer_inj_cells);
848  const ADB mob_others = subset(mob_perfcells[water_pos], other_well_cells);
849 
850  mob_polymer_inj = mob_polymer_inj / inv_wat_eff_visc / fully_mixing_visc;
851  mob_perfcells[water_pos] = superset(mob_polymer_inj, polymer_inj_cells, nperf) +
852  superset(mob_others, other_well_cells, nperf);
853  }
854  }
855  }
856 
857 } // namespace Opm
858 
859 #endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
void computeWaterShearVelocityFaces(const V &transi, const std::vector< ADB > &phasePressure, const SolutionState &state, std::vector< double > &water_vel, std::vector< double > &visc_mult)
Computing the water velocity without shear-thinning for the cell faces.
Definition: BlackoilPolymerModel_impl.hpp:599
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
Definition: PolymerPropsAd.hpp:32
SimulatorReport assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilPolymerModel_impl.hpp:496
int numPhases() const
Definition: BlackoilPropsAdFromDeck.cpp:225
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
Class for handling the standard well model.
Definition: StandardWells.hpp:51
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilPolymerModel_impl.hpp:409
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilPolymerModel_impl.hpp:124
void afterStep(const SimulatorTimerInterface &timer, ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilPolymerModel_impl.hpp:141
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
std::vector< ADB > select(const std::vector< ADB > &xc) const
Apply selector to multiple per-cell quantities.
Definition: AutoDiffHelpers.hpp:230
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
BlackoilPolymerModel(const typename Base::ModelParameters &param, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const PolymerPropsAd &polymer_props_ad, const StandardWells &well_model, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclipseState, const bool has_disgas, const bool has_vapoil, const bool has_polymer, const bool has_plyshlog, const bool has_shrate, const std::vector< double > &wells_rep_radius, const std::vector< double > &wells_perf_length, const std::vector< double > &wells_bore_diameter, const bool terminal_output)
Construct the model.
Definition: BlackoilPolymerModel_impl.hpp:77
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
A model implementation for three-phase black oil with polymer.
Definition: BlackoilPolymerModel.hpp:45
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
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
Definition: AutoDiffHelpers.hpp:623
Upwind selection in absence of counter-current flow (i.e., without effects of gravity and/or capillar...
Definition: AutoDiffHelpers.hpp:181
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
void computeWaterShearVelocityWells(const SolutionState &state, WellState &xw, const ADB &cq_sw, std::vector< double > &water_vel_wells, std::vector< double > &visc_mult_wells)
Computing the water velocity without shear-thinning for the well perforations based on the water flux...
Definition: BlackoilPolymerModel_impl.hpp:715