BlackoilReorderingTransportModel.hpp
1 /*
2  Copyright 2016 SINTEF ICT, Applied Mathematics.
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_BLACKOILREORDERINGTRANSPORTMODEL_HEADER_INCLUDED
21 #define OPM_BLACKOILREORDERINGTRANSPORTMODEL_HEADER_INCLUDED
22 
23 #include <opm/autodiff/BlackoilModelBase.hpp>
24 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
25 #include <opm/autodiff/BlackoilModelParameters.hpp>
26 #include <opm/autodiff/DebugTimeReport.hpp>
27 #include <opm/autodiff/multiPhaseUpwind.hpp>
28 #include <opm/core/grid.h>
29 #include <opm/core/transport/reorder/reordersequence.h>
30 #include <opm/core/simulator/BlackoilState.hpp>
31 
32 #include <opm/autodiff/BlackoilTransportModel.hpp>
33 
34 namespace Opm {
35 
36 
37  namespace detail
38  {
39 
40  template <typename Scalar>
42  {
43  Scalar operator()(double value, int index)
44  {
45  return Scalar::createVariable(value, index);
46  }
47  };
48 
49 
50 
51  template <>
52  struct CreateVariable<double>
53  {
54  double operator()(double value, int)
55  {
56  return value;
57  }
58  };
59 
60 
61 
62  template <typename Scalar>
64  {
65  Scalar operator()(double value)
66  {
67  return Scalar::createConstant(value);
68  }
69  };
70 
71 
72 
73  template <>
74  struct CreateConstant<double>
75  {
76  double operator()(double value)
77  {
78  return value;
79  }
80  };
81 
82 
83 
84 
85  struct Connection
86  {
87  Connection(const int ind, const double s) : index(ind), sign(s) {}
88  int index;
89  double sign;
90  };
91 
92 
93 
94 
95  class Connections;
96 
97 
98 
100  {
101  public:
102  explicit ConnectivityGraph(const HelperOps& ops)
103  : grad_(ops.grad)
104  , div_(ops.div)
105  {
106  grad_ia_ = grad_.outerIndexPtr();
107  grad_ja_ = grad_.innerIndexPtr();
108  grad_sign_ = grad_.valuePtr();
109  div_ia_ = div_.outerIndexPtr();
110  div_ja_ = div_.innerIndexPtr();
111  div_sign_ = div_.valuePtr();
112  }
113 
114  Connections cellConnections(const int cell) const;
115 
116  std::array<int, 2> connectionCells(const int connection) const
117  {
118  const int pos = div_ia_[connection];
119  assert(div_ia_[connection + 1] == pos + 2);
120  const double sign1 = div_sign_[pos];
121  assert(div_sign_[pos + 1] == -sign1);
122  if (sign1 > 0.0) {
123  return {{ div_ja_[pos], div_ja_[pos + 1] }};
124  } else {
125  return {{ div_ja_[pos + 1], div_ja_[pos] }};
126  }
127  }
128 
129  private:
130  friend class Connections;
131  typedef HelperOps::M M;
132  const M& grad_;
133  const M& div_;
134  const int* grad_ia_;
135  const int* grad_ja_;
136  const double* grad_sign_;
137  const int* div_ia_;
138  const int* div_ja_;
139  const double* div_sign_;
140  };
141 
142 
144  {
145  public:
146  Connections(const ConnectivityGraph& cg, const int cell) : cg_(cg), cell_(cell) {}
147  int size() const
148  {
149  return cg_.grad_ia_[cell_ + 1] - cg_.grad_ia_[cell_];
150  }
151  class Iterator
152  {
153  public:
154  Iterator(const Connections& c, const int index) : c_(c), index_(index) {}
155  Iterator& operator++()
156  {
157  ++index_;
158  return *this;
159  }
160  bool operator!=(const Iterator& other)
161  {
162  assert(&c_ == &other.c_);
163  return index_ != other.index_;
164  }
165  Connection operator*()
166  {
167  assert(index_ >= 0 && index_ < c_.size());
168  const int pos = c_.cg_.grad_ia_[c_.cell_] + index_;
169  return Connection(c_.cg_.grad_ja_[pos], -c_.cg_.grad_sign_[pos]); // Note the minus sign!
170  }
171  private:
172  const Connections& c_;
173  int index_;
174  };
175  Iterator begin() const { return Iterator(*this, 0); }
176  Iterator end() const { return Iterator(*this, size()); }
177  private:
178  friend class Iterator;
179  const ConnectivityGraph& cg_;
180  const int cell_;
181  };
182 
183 
184  inline Connections ConnectivityGraph::cellConnections(const int cell) const
185  {
186  return Connections(*this, cell);
187  }
188 
189 
190 
191  } // namespace detail
192 
193 
194 
195 
196 
197 
198 
200  template<class Grid, class WellModel>
202  : public BlackoilModelBase<Grid, WellModel, BlackoilReorderingTransportModel<Grid, WellModel> >
203  {
204  public:
206  friend Base;
207 
208  typedef typename Base::ReservoirState ReservoirState;
209  typedef typename Base::WellState WellState;
210  typedef typename Base::SolutionState SolutionState;
211  typedef typename Base::V V;
212 
213 
228  BlackoilReorderingTransportModel(const typename Base::ModelParameters& param,
229  const Grid& grid,
230  const BlackoilPropsAdFromDeck& fluid,
231  const DerivedGeology& geo,
232  const RockCompressibility* rock_comp_props,
233  const StandardWells& std_wells,
234  const NewtonIterationBlackoilInterface& linsolver,
235  std::shared_ptr<const EclipseState> eclState,
236  const bool has_disgas,
237  const bool has_vapoil,
238  const bool terminal_output)
239  : Base(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
240  eclState, has_disgas, has_vapoil, terminal_output)
241  , graph_(Base::ops_)
242  , props_(dynamic_cast<const BlackoilPropsAdFromDeck&>(fluid)) // TODO: remove the need for this cast.
243  , state0_{ ReservoirState(0, 0, 0), WellState(), V(), V() }
244  , state_{ ReservoirState(0, 0, 0), WellState(), V(), V() }
245  , tr_model_(param, grid, fluid, geo, rock_comp_props, std_wells, linsolver,
246  eclState, has_disgas, has_vapoil, terminal_output)
247  {
248  // Set up the common parts of the mass balance equations
249  // for each active phase.
250  const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
251  const V trans_nnc = ops_.nnc_trans;
252  trans_all_ = V::Zero(transi.size() + trans_nnc.size());
253  trans_all_ << transi, trans_nnc;
254  gdz_ = geo_.gravity()[2] * (ops_.grad * geo_.z().matrix());
255  rhos_ = DataBlock::Zero(ops_.div.rows(), 3);
256  rhos_.col(Water) = props_.surfaceDensity(Water, Base::cells_);
257  rhos_.col(Oil) = props_.surfaceDensity(Oil, Base::cells_);
258  rhos_.col(Gas) = props_.surfaceDensity(Gas, Base::cells_);
259  }
260 
261 
262 
263 
264 
265  void prepareStep(const SimulatorTimerInterface& timer,
266  const ReservoirState& reservoir_state,
267  const WellState& well_state)
268  {
269  tr_model_.prepareStep(timer, reservoir_state, well_state);
270  Base::prepareStep(timer, reservoir_state, well_state);
271  Base::param_.solve_welleq_initially_ = false;
272  state0_.reservoir_state = reservoir_state;
273  state0_.well_state = well_state;
274  // Since (reference) pressure is constant, porosity and transmissibility multipliers can
275  // be computed just once.
276  const std::vector<double>& p = reservoir_state.pressure();
277  state0_.tr_mult = Base::transMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
278  state0_.pv_mult = Base::poroMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
279  const int num_cells = p.size();
280  cstate0_.resize(num_cells);
281  for (int cell = 0; cell < num_cells; ++cell) {
282  computeCellState(cell, state0_, cstate0_[cell]);
283  }
284  cstate_ = cstate0_;
285  }
286 
287 
288 
289 
290 
291  template <class NonlinearSolverType>
292  SimulatorReport nonlinearIteration(const int iteration,
293  const SimulatorTimerInterface& timer,
294  NonlinearSolverType& nonlinear_solver,
295  ReservoirState& reservoir_state,
296  const WellState& well_state)
297  {
298  // Extract reservoir and well fluxes and state.
299  {
300  DebugTimeReport tr("Extracting fluxes");
301  extractFluxes(reservoir_state, well_state);
302  extractState(reservoir_state, well_state);
303  }
304 
305  // Compute cell ordering based on total flux.
306  {
307  DebugTimeReport tr("Topological sort");
308  computeOrdering();
309  }
310 
311  // Solve in every component (cell or block of cells), in order.
312  {
313  DebugTimeReport tr("Solving all components");
314  for (int ii = 0; ii < 5; ++ii) {
315  DebugTimeReport tr2("Solving components single sweep.");
316  solveComponents();
317  }
318  }
319 
320  // Update states for output.
321  reservoir_state = state_.reservoir_state;
322 
323  // Assemble with other model,
324  {
325  auto rs = reservoir_state;
326  auto ws = well_state;
327  tr_model_.nonlinearIteration(/*iteration*/ 0, timer, nonlinear_solver, rs, ws);
328  }
329 
330  // Create report and exit.
331  SimulatorReport report;
332  report.converged = true;
333  return report;
334  }
335 
336 
337 
338 
339 
340  void afterStep(const SimulatorTimerInterface& /* timer */,
341  const ReservoirState& /* reservoir_state */,
342  const WellState& /* well_state */)
343  {
344  // Does nothing in this model.
345  }
346 
347  using Base::numPhases;
348 
349 
350  protected:
351 
352  // ============ Types ============
353 
354  using Vec2 = Dune::FieldVector<double, 2>;
355  using Mat22 = Dune::FieldMatrix<double, 2, 2>;
356  using Eval = DenseAd::Evaluation<double, 2>;
357 
358 
359 
360 
361  struct State
362  {
363  ReservoirState reservoir_state;
364  WellState well_state;
365  V tr_mult;
366  V pv_mult;
367  };
368 
369 
370 
371 
372 
373  template <typename ScalarT>
374  struct CellState
375  {
376  using Scalar = ScalarT;
377 
378  Scalar s[3];
379  Scalar rs;
380  Scalar rv;
381  Scalar p[3];
382  Scalar kr[3];
383  Scalar pc[3];
384  Scalar temperature;
385  Scalar mu[3];
386  Scalar b[3];
387  Scalar lambda[3];
388  Scalar rho[3];
389  Scalar rssat;
390  Scalar rvsat;
391 
392  // Implement interface used for opm-material properties.
393  const Scalar& saturation(int phaseIdx) const
394  {
395  return s[phaseIdx];
396  }
397 
398  template <typename T>
399  CellState<T> flatten() const
400  {
401  return CellState<T>{
402  { s[0].value(), s[1].value(), s[2].value() },
403  rs.value(),
404  rv.value(),
405  { p[0].value(), p[1].value(), p[2].value() },
406  { kr[0].value(), kr[1].value(), kr[2].value() },
407  { pc[0].value(), pc[1].value(), pc[2].value() },
408  temperature.value(),
409  { mu[0].value(), mu[1].value(), mu[2].value() },
410  { b[0].value(), b[1].value(), b[2].value() },
411  { lambda[0].value(), lambda[1].value(), lambda[2].value() },
412  { rho[0].value(), rho[1].value(), rho[2].value() },
413  rssat.value(),
414  rvsat.value()
415  };
416  }
417  };
418 
419 
420 
421  // ============ Data members ============
422 
423  using Base::grid_;
424  using Base::geo_;
425  using Base::ops_;
426 
427  const detail::ConnectivityGraph graph_;
428 
429  const BlackoilPropsAdFromDeck& props_;
430 
431  State state0_;
432  State state_;
433 
434  std::vector<CellState<double>> cstate0_;
435  std::vector<CellState<double>> cstate_;
436 
437  V total_flux_;
438  V total_wellperf_flux_;
439  DataBlock comp_wellperf_flux_;
440  V total_wellflux_cell_;
441  V oil_wellflux_cell_;
442  V gas_wellflux_cell_;
443  std::vector<int> sequence_;
444  std::vector<int> components_;
445  V trans_all_;
446  V gdz_;
447  DataBlock rhos_;
448 
449  std::array<double, 2> max_abs_dx_;
450  std::array<int, 2> max_abs_dx_cell_;
451 
452  // TODO: remove this, for debug only.
454 
455 
456  // ============ Member functions ============
457 
458 
459 
460 
461 
462  template <typename Scalar>
463  void computeCellState(const int cell, const State& state, CellState<Scalar>& cstate) const
464  {
465  assert(numPhases() == 3); // I apologize for this to my future self, that will have to fix it.
466 
467  // Extract from state and props.
468  const auto hcstate = state.reservoir_state.hydroCarbonState()[cell];
469  const bool is_sg = (hcstate == HydroCarbonState::GasAndOil);
470  const bool is_rs = (hcstate == HydroCarbonState::OilOnly);
471  const bool is_rv = (hcstate == HydroCarbonState::GasOnly);
472  const double swval = state.reservoir_state.saturation()[3*cell + Water];
473  const double sgval = state.reservoir_state.saturation()[3*cell + Gas];
474  const double rsval = state.reservoir_state.gasoilratio()[cell];
475  const double rvval = state.reservoir_state.rv()[cell];
476  const double poval = state.reservoir_state.pressure()[cell];
477  const int pvt_region = props_.pvtRegions()[cell];
478 
479  // Property functions.
480  const auto& waterpvt = props_.waterProps();
481  const auto& oilpvt = props_.oilProps();
482  const auto& gaspvt = props_.gasProps();
483  const auto& satfunc = props_.materialLaws();
484 
485  // Create saturation and composition variables.
488  cstate.s[Water] = variable(swval, 0);
489  cstate.s[Gas] = is_sg ? variable(sgval, 1) : constant(sgval);
490  cstate.s[Oil] = 1.0 - cstate.s[Water] - cstate.s[Gas];
491  cstate.rs = is_rs ? variable(rsval, 1) : constant(rsval);
492  cstate.rv = is_rv ? variable(rvval, 1) : constant(rvval);
493 
494  // Compute relative permeabilities amd capillary pressures.
495  const auto& params = satfunc.materialLawParams(cell);
496  typedef BlackoilPropsAdFromDeck::MaterialLawManager::MaterialLaw MaterialLaw;
497  MaterialLaw::relativePermeabilities(cstate.kr, params, cstate);
498  MaterialLaw::capillaryPressures(cstate.pc, params, cstate);
499 
500  // Compute phase pressures.
501  cstate.p[Oil] = constant(poval);
502  cstate.p[Water] = cstate.p[Oil] + cstate.pc[Water]; // pcow = pw - po (!) [different from old convention]
503  cstate.p[Gas] = cstate.p[Oil] + cstate.pc[Gas]; // pcog = pg - po
504 
505  // Compute PVT properties.
506  cstate.temperature = constant(0.0); // Temperature is not used.
507  cstate.mu[Water] = waterpvt.viscosity(pvt_region, cstate.temperature, cstate.p[Water]);
508  cstate.mu[Oil] = is_sg
509  ? oilpvt.saturatedViscosity(pvt_region, cstate.temperature, cstate.p[Oil])
510  : oilpvt.viscosity(pvt_region, cstate.temperature, cstate.p[Oil], cstate.rs);
511  cstate.mu[Gas] = is_sg
512  ? gaspvt.saturatedViscosity(pvt_region, cstate.temperature, cstate.p[Gas])
513  : gaspvt.viscosity(pvt_region, cstate.temperature, cstate.p[Gas], cstate.rv);
514  cstate.b[Water] = waterpvt.inverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Water]);
515  cstate.b[Oil] = is_sg
516  ? oilpvt.saturatedInverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Oil])
517  : oilpvt.inverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Oil], cstate.rs);
518  cstate.b[Gas] = is_sg
519  ? gaspvt.saturatedInverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Gas])
520  : gaspvt.inverseFormationVolumeFactor(pvt_region, cstate.temperature, cstate.p[Gas], cstate.rv);
521 
522  // Compute mobilities.
523  for (int phase = 0; phase < 3; ++phase) {
524  cstate.lambda[phase] = cstate.kr[phase] / cstate.mu[phase];
525  }
526 
527  // Compute densities.
528  cstate.rho[Water] = rhos_(cell, Water) * cstate.b[Water];
529  cstate.rho[Oil] = (rhos_(cell, Oil) + cstate.rs*rhos_(cell, Gas)) * cstate.b[Oil]; // TODO: check that this is correct
530  cstate.rho[Gas] = (rhos_(cell, Gas) + cstate.rv*rhos_(cell, Oil)) * cstate.b[Gas];
531 
532  // Compute saturated rs and rv factors.
533  cstate.rssat = oilpvt.saturatedGasDissolutionFactor(pvt_region, cstate.temperature, cstate.p[Oil]);
534  cstate.rvsat = gaspvt.saturatedOilVaporizationFactor(pvt_region, cstate.temperature, cstate.p[Gas]);
535  // TODO: add vaporization controls such as in BlackoilPropsAdFromDeck::applyVap().
536  }
537 
538 
539 
540 
541  void extractFluxes(const ReservoirState& reservoir_state,
542  const WellState& well_state)
543  {
544  // Input face fluxes are for interior faces + nncs.
545  total_flux_ = Eigen::Map<const V>(reservoir_state.faceflux().data(),
546  reservoir_state.faceflux().size());
547  total_wellperf_flux_ = Eigen::Map<const V>(well_state.perfRates().data(),
548  well_state.perfRates().size());
549  comp_wellperf_flux_ = Eigen::Map<const DataBlock>(well_state.perfPhaseRates().data(),
550  well_state.perfRates().size(),
551  numPhases());
552  const int num_cells = reservoir_state.pressure().size();
553  total_wellflux_cell_ = superset(total_wellperf_flux_, Base::wellModel().wellOps().well_cells, num_cells);
554  assert(Base::numPhases() == 3);
555  V oilflux = comp_wellperf_flux_.col(1);
556  V gasflux = comp_wellperf_flux_.col(2);
557  oil_wellflux_cell_ = superset(oilflux, Base::wellModel().wellOps().well_cells, num_cells);
558  gas_wellflux_cell_ = superset(gasflux, Base::wellModel().wellOps().well_cells, num_cells);
559  assert(numPhases() * well_state.perfRates().size() == well_state.perfPhaseRates().size());
560  }
561 
562 
563 
564 
565 
566  void extractState(const ReservoirState& reservoir_state,
567  const WellState& well_state)
568  {
569  state_.reservoir_state = reservoir_state;
570  state_.well_state = well_state;
571  const std::vector<double>& p = reservoir_state.pressure();
572  state_.tr_mult = Base::transMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
573  state_.pv_mult = Base::poroMult(ADB::constant(Eigen::Map<const V>(p.data(), p.size()))).value();
574  }
575 
576 
577 
578 
579 
580  void computeOrdering()
581  {
582  assert(!geo_.nnc().hasNNC()); // TODO: support compute_sequence() with grid + nnc.
583  static_assert(std::is_same<Grid, UnstructuredGrid>::value,
584  "compute_sequence() is written in C and therefore requires an UnstructuredGrid, "
585  "it must be rewritten to use other grid classes such as CpGrid");
586  using namespace Opm::AutoDiffGrid;
587  const int num_cells = numCells(grid_);
588  sequence_.resize(num_cells);
589  components_.resize(num_cells + 1); // max possible size
590  int num_components = -1;
591 
592 
593  using namespace Opm::AutoDiffGrid;
594  const int num_faces = numFaces(grid_);
595  V flux_on_all_faces = superset(total_flux_, ops_.internal_faces, num_faces);
596  compute_sequence(&grid_, flux_on_all_faces.data(), sequence_.data(), components_.data(), &num_components);
597  OpmLog::debug(std::string("Number of components: ") + std::to_string(num_components));
598  components_.resize(num_components + 1); // resize to fit actually used part
599  }
600 
601 
602 
603 
604  void solveComponents()
605  {
606  // Zero the max changed.
607  max_abs_dx_[0] = 0.0;
608  max_abs_dx_[1] = 0.0;
609  max_abs_dx_cell_[0] = -1;
610  max_abs_dx_cell_[1] = -1;
611 
612  // Solve the equations.
613  const int num_components = components_.size() - 1;
614  for (int comp = 0; comp < num_components; ++comp) {
615  const int comp_size = components_[comp + 1] - components_[comp];
616  if (comp_size == 1) {
617  solveSingleCell(sequence_[components_[comp]]);
618  } else {
619  solveMultiCell(comp_size, &sequence_[components_[comp]]);
620  }
621  }
622 
623  // Log the max change.
624  {
625  std::ostringstream os;
626  os << "=== Max abs dx[0]: " << max_abs_dx_[0] << " (cell " << max_abs_dx_cell_[0]
627  <<") dx[1]: " << max_abs_dx_[1] << " (cell " << max_abs_dx_cell_[1] << ")";
628  OpmLog::debug(os.str());
629  }
630  }
631 
632 
633 
634 
635 
636  void solveSingleCell(const int cell)
637  {
638 
639  Vec2 res;
640  Mat22 jac;
641  assembleSingleCell(cell, res, jac);
642 
643  // Newton loop.
644  int iter = 0;
645  const int max_iter = 100;
646  double relaxation = 1.0;
647  while (!getConvergence(cell, res) && iter < max_iter) {
648  Vec2 dx;
649  jac.solve(dx, res);
650  dx *= relaxation;
651  // const auto hcstate_old = state_.reservoir_state.hydroCarbonState()[cell];
652  updateState(cell, -dx);
653  // const auto hcstate = state_.reservoir_state.hydroCarbonState()[cell];
654  assembleSingleCell(cell, res, jac);
655  ++iter;
656  if (iter > 10) {
657  relaxation = 0.85;
658  if (iter > 15) {
659  relaxation = 0.70;
660  }
661  if (iter > 20) {
662  relaxation = 0.55;
663  }
664  if (iter > 25) {
665  relaxation = 0.40;
666  }
667  if (iter > 30) {
668  relaxation = 0.25;
669  }
670  // std::ostringstream os;
671  // os << "Iteration " << iter << " in cell " << cell << ", residual = " << res
672  // << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas]
673  // << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << "}, dx = " << dx << ", hcstate: " << hcstate_old << " -> " << hcstate;
674  // OpmLog::debug(os.str());
675  }
676  }
677  if (iter == max_iter) {
678  std::ostringstream os;
679  os << "Failed to converge in cell " << cell << ", residual = " << res
680  << ", cell values { s = ( " << cstate_[cell].s[Water] << ", " << cstate_[cell].s[Oil] << ", " << cstate_[cell].s[Gas]
681  << " ), rs = " << cstate_[cell].rs << ", rv = " << cstate_[cell].rv << " }";
682  OpmLog::debug(os.str());
683  }
684  }
685 
686 
687 
688 
689 
690  void solveMultiCell(const int comp_size, const int* cell_array)
691  {
692  // OpmLog::warning("solveMultiCell", "solveMultiCell() called with component size " + std::to_string(comp_size));
693  for (int ii = 0; ii < comp_size; ++ii) {
694  solveSingleCell(cell_array[ii]);
695  }
696  }
697 
698 
699 
700 
701  template <typename Scalar>
702  Scalar oilAccumulation(const CellState<Scalar>& cs)
703  {
704  return cs.b[Oil]*cs.s[Oil] + cs.rv*cs.b[Gas]*cs.s[Gas];
705  }
706 
707 
708 
709 
710  template <typename Scalar>
711  Scalar gasAccumulation(const CellState<Scalar>& cs)
712  {
713  return cs.b[Gas]*cs.s[Gas] + cs.rs*cs.b[Oil]*cs.s[Oil];
714  }
715 
716 
717 
718 
719  void applyThresholdPressure(const int connection, Eval& dp)
720  {
721  const double thres_press = Base::threshold_pressures_by_connection_[connection];
722  if (std::fabs(dp.value()) < thres_press) {
723  dp.setValue(0.0);
724  } else {
725  dp -= dp.value() > 0.0 ? thres_press : -thres_press;
726  }
727  }
728 
729 
730 
731 
732  void assembleSingleCell(const int cell, Vec2& res, Mat22& jac)
733  {
734  assert(numPhases() == 3); // I apologize for this to my future self, that will have to fix it.
735 
736  CellState<Eval> st;
737  computeCellState(cell, state_, st);
738  cstate_[cell] = st.template flatten<double>();
739 
740  // Accumulation terms.
741  const double pvm0 = state0_.pv_mult[cell];
742  const double pvm = state_.pv_mult[cell];
743  const double ao0 = oilAccumulation(cstate0_[cell]) * pvm0;
744  const Eval ao = oilAccumulation(st) * pvm;
745  const double ag0 = gasAccumulation(cstate0_[cell]) * pvm0;
746  const Eval ag = gasAccumulation(st) * pvm;
747 
748  // Flux terms.
749  Eval div_oilflux = Eval::createConstant(0.0);
750  Eval div_gasflux = Eval::createConstant(0.0);
751  for (auto conn : graph_.cellConnections(cell)) {
752  auto conn_cells = graph_.connectionCells(conn.index);
753  const int from = conn_cells[0];
754  const int to = conn_cells[1];
755  if (from < 0 || to < 0) {
756  continue; // Boundary.
757  }
758  assert((from == cell) == (conn.sign > 0.0));
759  const int other = from == cell ? to : from;
760  const double vt = conn.sign * total_flux_[conn.index];
761  const double gdz = conn.sign * gdz_[conn.index];
762 
763  // From this point, we treat everything about this
764  // connection as going from 'cell' to 'other'. Since
765  // we don't want derivatives from the 'other' cell to
766  // participate in the solution, we use the constant
767  // values from cstate_[other].
768  Eval dh[3];
769  Eval dh_sat[3];
770  const Eval grad_oil_press = cstate_[other].p[Oil] - st.p[Oil];
771  for (int phase : { Water, Oil, Gas }) {
772  const Eval gradp = cstate_[other].p[phase] - st.p[phase];
773  const Eval rhoavg = 0.5 * (st.rho[phase] + cstate_[other].rho[phase]);
774  dh[phase] = gradp - rhoavg * gdz;
775  if (Base::use_threshold_pressure_) {
776  applyThresholdPressure(conn.index, dh[phase]);
777  }
778  dh_sat[phase] = grad_oil_press - dh[phase];
779  }
780  const double tran = trans_all_[conn.index]; // TODO: include tr_mult effect.
781  const auto& m1 = st.lambda;
782  const auto& m2 = cstate_[other].lambda;
783  const auto upw = connectionMultiPhaseUpwind({{ dh_sat[Water].value(), dh_sat[Oil].value(), dh_sat[Gas].value() }},
784  {{ m1[Water].value(), m1[Oil].value(), m1[Gas].value() }},
785  {{ m2[Water], m2[Oil], m2[Gas] }},
786  tran, vt);
787  // if (upw[0] != upw[1] || upw[1] != upw[2]) {
788  // OpmLog::debug("Detected countercurrent flow between cells " + std::to_string(from) + " and " + std::to_string(to));
789  // }
790  Eval b[3];
791  Eval mob[3];
792  Eval tot_mob = Eval::createConstant(0.0);
793  for (int phase : { Water, Oil, Gas }) {
794  b[phase] = upw[phase] > 0.0 ? st.b[phase] : cstate_[other].b[phase];
795  mob[phase] = upw[phase] > 0.0 ? m1[phase] : m2[phase];
796  tot_mob += mob[phase];
797  }
798  Eval rs = upw[Oil] > 0.0 ? st.rs : cstate_[other].rs;
799  Eval rv = upw[Gas] > 0.0 ? st.rv : cstate_[other].rv;
800 
801  Eval flux[3];
802  for (int phase : { Oil, Gas }) {
803  Eval gflux = Eval::createConstant(0.0);
804  for (int other_phase : { Water, Oil, Gas }) {
805  if (phase != other_phase) {
806  gflux += mob[other_phase] * (dh_sat[phase] - dh_sat[other_phase]);
807  }
808  }
809  flux[phase] = b[phase] * (mob[phase] / tot_mob) * (vt + tran*gflux);
810  }
811  div_oilflux += flux[Oil] + rv*flux[Gas];
812  div_gasflux += flux[Gas] + rs*flux[Oil];
813  }
814 
815  // Well fluxes.
816  if (total_wellflux_cell_[cell] > 0.0) {
817  // Injecting perforation. Use given phase rates.
818  assert(oil_wellflux_cell_[cell] >= 0.0);
819  assert(gas_wellflux_cell_[cell] >= 0.0);
820  div_oilflux -= oil_wellflux_cell_[cell];
821  div_gasflux -= gas_wellflux_cell_[cell];
822  } else if (total_wellflux_cell_[cell] < 0.0) {
823  // Producing perforation. Use total rate and fractional flow.
824  Eval totmob = st.lambda[Water] + st.lambda[Oil] + st.lambda[Gas];
825  Eval oilflux = st.b[Oil] * (st.lambda[Oil]/totmob) * total_wellflux_cell_[cell];
826  Eval gasflux = st.b[Gas] * (st.lambda[Gas]/totmob) * total_wellflux_cell_[cell];
827  div_oilflux -= (oilflux + st.rv * gasflux);
828  div_gasflux -= (gasflux + st.rs * oilflux);
829  }
830 
831  const Eval oileq = Base::pvdt_[cell]*(ao - ao0) + div_oilflux;
832  const Eval gaseq = Base::pvdt_[cell]*(ag - ag0) + div_gasflux;
833 
834  res[0] = oileq.value();
835  res[1] = gaseq.value();
836  jac[0][0] = oileq.derivative(0);
837  jac[0][1] = oileq.derivative(1);
838  jac[1][0] = gaseq.derivative(0);
839  jac[1][1] = gaseq.derivative(1);
840  }
841 
842 
843 
844 
845 
846  bool getConvergence(const int cell, const Vec2& res)
847  {
848  const double tol = 1e-7;
849  // Compute scaled residuals (scaled like saturations).
850  double sres[] = { res[0] / (cstate_[cell].b[Oil] * Base::pvdt_[cell]),
851  res[1] / (cstate_[cell].b[Gas] * Base::pvdt_[cell]) };
852  return std::fabs(sres[0]) < tol && std::fabs(sres[1]) < tol;
853  }
854 
855 
856 
857 
858  void updateState(const int cell,
859  const Vec2& dx)
860  {
861  if (std::fabs(dx[0]) > max_abs_dx_[0]) {
862  max_abs_dx_cell_[0] = cell;
863  }
864  if (std::fabs(dx[1]) > max_abs_dx_[1]) {
865  max_abs_dx_cell_[1] = cell;
866  }
867  max_abs_dx_[0] = std::max(max_abs_dx_[0], std::fabs(dx[0]));
868  max_abs_dx_[1] = std::max(max_abs_dx_[1], std::fabs(dx[1]));
869 
870  // Get saturation updates.
871  const double dsw = dx[0];
872  double dsg = 0.0;
873  auto& hcstate = state_.reservoir_state.hydroCarbonState()[cell];
874  if (hcstate == HydroCarbonState::GasAndOil) {
875  dsg = dx[1];
876  } else if (hcstate == HydroCarbonState::GasOnly) {
877  dsg = -dsw;
878  }
879  const double dso = -(dsw + dsg);
880 
881  // Handle too large saturation changes.
882  const double maxval = std::max(std::fabs(dsw), std::max(std::fabs(dso), std::fabs(dsg)));
883  const double sfactor = std::min(1.0, Base::dsMax() / maxval);
884  double* s = state_.reservoir_state.saturation().data() + 3*cell;
885  s[Water] += sfactor*dsw;
886  s[Gas] += sfactor*dsg;
887  s[Oil] = 1.0 - s[Water] - s[Gas];
888 
889  // Handle < 0 saturations.
890  for (int phase : { Gas, Oil, Water }) { // TODO: check if ordering here is significant
891  if (s[phase] < 0.0) {
892  for (int other_phase : { Water, Oil, Gas }) {
893  if (phase != other_phase) {
894  s[other_phase] /= (1.0 - s[phase]);
895  }
896  }
897  s[phase] = 0.0;
898  }
899  }
900 
901  // Update rs.
902  double& rs = state_.reservoir_state.gasoilratio()[cell];
903  const double rs_old = rs;
904  if (hcstate == HydroCarbonState::OilOnly) {
905  // const double max_allowed_change = std::fabs(rs_old) * Base::drMaxRel();
906  const double drs = dx[1];
907  // const double factor = std::min(1.0, max_allowed_change / std::fabs(drs));
908  // rs += factor*drs;
909  rs += drs;
910  rs = std::max(rs, 0.0);
911  }
912 
913  // Update rv.
914  double& rv = state_.reservoir_state.rv()[cell];
915  const double rv_old = rv;
916  if (hcstate == HydroCarbonState::GasOnly) {
917  // const double max_allowed_change = std::fabs(rv_old) * Base::drMaxRel();
918  const double drv = dx[1];
919  // const double factor = std::min(1.0, max_allowed_change / std::fabs(drv));
920  // rv += factor*drv;
921  rv += drv;
922  rv = std::max(rv, 0.0);
923  }
924 
925  const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
926  const bool water_only = s[Water] > (1 - epsilon);
927  const auto old_hcstate = hcstate;
928  hcstate = HydroCarbonState::GasAndOil;
929  // sg <-> rs transition.
930  {
931  const double rssat_old = cstate_[cell].rssat;
932  const double rssat = rssat_old; // TODO: This is no longer true with vaporization controls
933  const bool is_rs = old_hcstate == HydroCarbonState::OilOnly;
934  const bool has_gas = (s[Gas] > 0.0 && !is_rs);
935  const bool gas_vaporized = ( (rs > rssat * (1+epsilon) && is_rs ) && (rs_old > rssat_old * (1-epsilon)) );
936  if (water_only || has_gas || gas_vaporized) {
937  rs = rssat;
938  } else {
939  hcstate = HydroCarbonState::OilOnly;
940  }
941  }
942 
943  // sg <-> rv transition.
944  {
945  const double rvsat_old = cstate_[cell].rvsat;
946  const double rvsat = rvsat_old; // TODO: This is no longer true with vaporization controls
947  const bool is_rv = old_hcstate == HydroCarbonState::GasOnly;
948  const bool has_oil = (s[Oil] > 0.0 && !is_rv);
949  const bool oil_condensed = ( (rv > rvsat * (1+epsilon) && is_rv) && (rv_old > rvsat_old * (1-epsilon)) );
950  if (water_only || has_oil || oil_condensed) {
951  rv = rvsat;
952  } else {
953  hcstate = HydroCarbonState::GasOnly;
954  }
955  }
956  }
957  };
958 
959 
960 
961 
962 
963 
964 
965 
967  template <class Grid, class WellModel>
969  {
970  typedef BlackoilState ReservoirState;
974  };
975 
976 } // namespace Opm
977 
978 
979 
980 
981 #endif // OPM_BLACKOILREORDERINGTRANSPORTMODEL_HEADER_INCLUDED
Contains vectors and sparse matrices that represent subsets or operations on (AD or regular) vectors ...
Definition: AutoDiffHelpers.hpp:44
const MaterialLawManager & materialLaws() const
Direct access to lower-level saturation functions.
Definition: BlackoilPropsAdFromDeck.hpp:449
V surfaceDensity(const int phaseIdx, const Cells &cells) const
Densities of stock components at surface conditions.
Definition: BlackoilPropsAdFromDeck.cpp:242
V nnc_trans
The NNC transmissibilities.
Definition: AutoDiffHelpers.hpp:72
M div
Extract for each cell the sum of its adjacent interior faces&#39; (signed) values.
Definition: AutoDiffHelpers.hpp:60
Definition: BlackoilReorderingTransportModel.hpp:151
BlackoilReorderingTransportModel(const typename Base::ModelParameters &param, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const StandardWells &std_wells, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output)
Construct the model.
Definition: BlackoilReorderingTransportModel.hpp:228
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
const OilPvt & oilProps() const
Direct access to lower-level oil pvt props.
Definition: BlackoilPropsAdFromDeck.hpp:437
Definition: BlackoilReorderingTransportModel.hpp:63
A model implementation for the transport equation in three-phase black oil.
Definition: BlackoilReorderingTransportModel.hpp:201
const GasPvt & gasProps() const
Direct access to lower-level gas pvt props.
Definition: BlackoilPropsAdFromDeck.hpp:443
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
std::array< double, 3 > connectionMultiPhaseUpwind(const std::array< double, 3 > &head_diff, const std::array< double, 3 > &mob1, const std::array< double, 3 > &mob2, const double transmissibility, const double flux)
Compute upwind directions for three-phase flow across a connection.
Definition: multiPhaseUpwind.cpp:31
Definition: BlackoilReorderingTransportModel.hpp:41
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
Class for handling the standard well model.
Definition: StandardWells.hpp:51
Definition: BlackoilReorderingTransportModel.hpp:143
Definition: BlackoilReorderingTransportModel.hpp:361
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilModelBase_impl.hpp:220
Definition: BlackoilReorderingTransportModel.hpp:374
Definition: GridHelpers.cpp:27
Traits to encapsulate the types used by classes using or extending this model.
Definition: BlackoilModelBase.hpp:60
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: BlackoilReorderingTransportModel.hpp:99
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
Struct for containing iteration variables.
Definition: DefaultBlackoilSolutionState.hpp:29
int rows() const
Returns number of rows in the matrix.
Definition: AutoDiffMatrix.hpp:565
Definition: BlackoilReorderingTransportModel.hpp:85
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelBase_impl.hpp:383
M grad
Extract for each face the difference of its adjacent cells&#39; values (second - first).
Definition: AutoDiffHelpers.hpp:56
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
const WaterPvt & waterProps() const
Direct access to lower-level water pvt props.
Definition: BlackoilPropsAdFromDeck.hpp:431
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
A model implementation for the transport equation in three-phase black oil.
Definition: BlackoilTransportModel.hpp:35