All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
RateConverter.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014, 2015 Statoil ASA.
4 
5  This file is part of the Open Porous Media Project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED
22 #define OPM_RATECONVERTER_HPP_HEADER_INCLUDED
23 
24 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
25 
26 #include <opm/core/props/BlackoilPhases.hpp>
27 #include <opm/core/simulator/BlackoilState.hpp>
28 #include <opm/core/utility/RegionMapping.hpp>
29 #include <opm/core/linalg/ParallelIstlInformation.hpp>
30 
31 #include <algorithm>
32 #include <cmath>
33 #include <memory>
34 #include <stdexcept>
35 #include <type_traits>
36 #include <unordered_map>
37 #include <utility>
38 #include <vector>
49 namespace Opm {
50  namespace RateConverter {
55  namespace Details {
56  namespace Select {
57  template <class RegionID, bool>
59  {
60  using type =
61  typename std::remove_reference<RegionID>::type &;
62  };
63 
64  template <class RegionID>
65  struct RegionIDParameter<RegionID, true>
66  {
67  using type = RegionID;
68  };
69  } // Select
70 
77  template<bool is_parallel>
79  {
91  std::tuple<double, double, double, double, int>
92  operator()(const std::vector<double>& pressure,
93  const std::vector<double>& temperature,
94  const std::vector<double>& rs,
95  const std::vector<double>& rv,
96  const std::vector<double>& ownership,
97  std::size_t cell){
98  if ( ownership[cell] )
99  {
100  return std::make_tuple(pressure[cell],
101  temperature[cell],
102  rs[cell],
103  rv[cell],
104  1);
105  }
106  else
107  {
108  return std::make_tuple(0, 0, 0, 0, 0);
109  }
110  }
111  };
112  template<>
114  {
115  std::tuple<double, double, double, double, int>
116  operator()(const std::vector<double>& pressure,
117  const std::vector<double>& temperature,
118  const std::vector<double>& rs,
119  const std::vector<double>& rv,
120  const std::vector<double>&,
121  std::size_t cell){
122  return std::make_tuple(pressure[cell],
123  temperature[cell],
124  rs[cell],
125  rv[cell],
126  1);
127  }
128  };
141  template <typename RegionId, class Attributes>
143  {
144  public:
149  using RegionID =
151  <RegionId, std::is_integral<RegionId>::value>::type;
152 
166  template <class RMap>
167  RegionAttributes(const RMap& rmap,
168  const Attributes& attr)
169  {
170  using VT = typename AttributeMap::value_type;
171 
172  for (const auto& r : rmap.activeRegions()) {
173  auto v = std::unique_ptr<Value>(new Value(attr));
174 
175  const auto stat = attr_.insert(VT(r, std::move(v)));
176 
177  if (stat.second) {
178  // New value inserted.
179  const auto& cells = rmap.cells(r);
180 
181  assert (! cells.empty());
182 
183  // Region's representative cell.
184  stat.first->second->cell_ = cells[0];
185  }
186  }
187  }
188 
196  int cell(const RegionID reg) const
197  {
198  return this->find(reg).cell_;
199  }
200 
209  const Attributes& attributes(const RegionID reg) const
210  {
211  return this->find(reg).attr_;
212  }
213 
222  Attributes& attributes(const RegionID reg)
223  {
224  return this->find(reg).attr_;
225  }
226 
227  private:
232  struct Value {
233  Value(const Attributes& attr)
234  : attr_(attr)
235  , cell_(-1)
236  {}
237 
238  Attributes attr_;
239  int cell_;
240  };
241 
242  using ID =
243  typename std::remove_reference<RegionId>::type;
244 
245  using AttributeMap =
246  std::unordered_map<ID, std::unique_ptr<Value>>;
247 
248  AttributeMap attr_;
249 
253  const Value& find(const RegionID reg) const
254  {
255  const auto& i = attr_.find(reg);
256 
257  if (i == attr_.end()) {
258  throw std::invalid_argument("Unknown region ID");
259  }
260 
261  return *i->second;
262  }
263 
267  Value& find(const RegionID reg)
268  {
269  const auto& i = attr_.find(reg);
270 
271  if (i == attr_.end()) {
272  throw std::invalid_argument("Unknown region ID");
273  }
274 
275  return *i->second;
276  }
277  };
278 
283  namespace PhaseUsed {
291  inline bool
292  water(const PhaseUsage& pu)
293  {
294  return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
295  }
296 
304  inline bool
305  oil(const PhaseUsage& pu)
306  {
307  return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
308  }
309 
317  inline bool
318  gas(const PhaseUsage& pu)
319  {
320  return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
321  }
322  } // namespace PhaseUsed
323 
328  namespace PhasePos {
337  inline int
338  water(const PhaseUsage& pu)
339  {
340  int p = -1;
341 
342  if (PhaseUsed::water(pu)) {
343  p = pu.phase_pos[ BlackoilPhases::Aqua ];
344  }
345 
346  return p;
347  }
348 
357  inline int
358  oil(const PhaseUsage& pu)
359  {
360  int p = -1;
361 
362  if (PhaseUsed::oil(pu)) {
363  p = pu.phase_pos[ BlackoilPhases::Liquid ];
364  }
365 
366  return p;
367  }
368 
377  inline int
378  gas(const PhaseUsage& pu)
379  {
380  int p = -1;
381 
382  if (PhaseUsed::gas(pu)) {
383  p = pu.phase_pos[ BlackoilPhases::Vapour ];
384  }
385 
386  return p;
387  }
388  } // namespace PhasePos
389  } // namespace Details
390 
406  template <class FluidSystem, class Region>
408  public:
416  SurfaceToReservoirVoidage(const PhaseUsage& phaseUsage,
417  const Region& region)
418  : phaseUsage_(phaseUsage)
419  , rmap_ (region)
420  , attr_ (rmap_, Attributes())
421  {
422  }
423 
439  void
440  defineState(const BlackoilState& state,
441  const boost::any& info = boost::any())
442  {
443 #if HAVE_MPI
444  if( info.type() == typeid(ParallelISTLInformation) )
445  {
446  const auto& ownership =
447  boost::any_cast<const ParallelISTLInformation&>(info)
448  .updateOwnerMask(state.pressure());
449  calcAverages<true>(state, info, ownership);
450  }
451  else
452 #endif
453  {
454  std::vector<double> dummyOwnership; // not actually used
455  calcAverages<false>(state, info, dummyOwnership);
456  }
457  }
458 
467  template <typename ElementContext, class EbosSimulator>
468  void defineState(const EbosSimulator& simulator)
469  {
470 
471  // create map from cell to region
472  // and set all attributes to zero
473  const auto& grid = simulator.gridManager().grid();
474  const unsigned numCells = grid.size(/*codim=*/0);
475  std::vector<int> cell2region(numCells, -1);
476  for (const auto& reg : rmap_.activeRegions()) {
477  for (const auto& cell : rmap_.cells(reg)) {
478  cell2region[cell] = reg;
479  }
480  auto& ra = attr_.attributes(reg);
481  ra.pressure = 0.0;
482  ra.temperature = 0.0;
483  ra.rs = 0.0;
484  ra.rv = 0.0;
485  ra.pv = 0.0;
486 
487  }
488 
489  ElementContext elemCtx( simulator );
490  const auto& gridView = simulator.gridView();
491  const auto& comm = gridView.comm();
492 
493  const auto& elemEndIt = gridView.template end</*codim=*/0>();
494  for (auto elemIt = gridView.template begin</*codim=*/0>();
495  elemIt != elemEndIt;
496  ++elemIt)
497  {
498 
499  const auto& elem = *elemIt;
500  if (elem.partitionType() != Dune::InteriorEntity)
501  continue;
502 
503  elemCtx.updatePrimaryStencil(elem);
504  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
505  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
506  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
507  const auto& fs = intQuants.fluidState();
508  // use pore volume weighted averages.
509  const double pv_cell =
510  simulator.model().dofTotalVolume(cellIdx)
511  * intQuants.porosity().value();
512 
513  // only count oil and gas filled parts of the domain
514  double hydrocarbon = 1.0;
515  const auto& pu = phaseUsage_;
516  if (Details::PhaseUsed::water(pu)) {
517  hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
518  }
519 
520  int reg = cell2region[cellIdx];
521  assert(reg >= 0);
522  auto& ra = attr_.attributes(reg);
523  auto& p = ra.pressure;
524  auto& T = ra.temperature;
525  auto& rs = ra.rs;
526  auto& rv = ra.rv;
527  auto& pv = ra.pv;
528 
529  // sum p, rs, rv, and T.
530  double hydrocarbonPV = pv_cell*hydrocarbon;
531  pv += hydrocarbonPV;
532  p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
533  rs += fs.Rs().value()*hydrocarbonPV;
534  rv += fs.Rv().value()*hydrocarbonPV;
535  T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
536  }
537 
538  for (const auto& reg : rmap_.activeRegions()) {
539  auto& ra = attr_.attributes(reg);
540  auto& p = ra.pressure;
541  auto& T = ra.temperature;
542  auto& rs = ra.rs;
543  auto& rv = ra.rv;
544  auto& pv = ra.pv;
545  // communicate sums
546  p = comm.sum(p);
547  T = comm.sum(T);
548  rs = comm.sum(rs);
549  rv = comm.sum(rv);
550  pv = comm.sum(pv);
551  // compute average
552  p /= pv;
553  T /= pv;
554  rs /= pv;
555  rv /= pv;
556  }
557  }
558 
564  typedef typename RegionMapping<Region>::RegionId RegionId;
565 
587  template <class Coeff>
588  void
589  calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
590  {
591  const auto& pu = phaseUsage_;
592  const auto& ra = attr_.attributes(r);
593  const double p = ra.pressure;
594  const double T = ra.temperature;
595 
596  const int iw = Details::PhasePos::water(pu);
597  const int io = Details::PhasePos::oil (pu);
598  const int ig = Details::PhasePos::gas (pu);
599 
600  std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
601 
602  if (Details::PhaseUsed::water(pu)) {
603  // q[w]_r = q[w]_s / bw
604 
605  const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
606 
607  coeff[iw] = 1.0 / bw;
608  }
609 
610  // Determinant of 'R' matrix
611  const double detR = 1.0 - (ra.rs * ra.rv);
612 
613  if (Details::PhaseUsed::oil(pu)) {
614  // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
615 
616  const double Rs = ra.rs;
617  const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
618  const double den = bo * detR;
619 
620  coeff[io] += 1.0 / den;
621 
622  if (Details::PhaseUsed::gas(pu)) {
623  coeff[ig] -= ra.rv / den;
624  }
625  }
626 
627  if (Details::PhaseUsed::gas(pu)) {
628  // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
629 
630  const double Rv = ra.rv;
631  const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
632  const double den = bg * detR;
633 
634  coeff[ig] += 1.0 / den;
635 
636  if (Details::PhaseUsed::oil(pu)) {
637  coeff[io] -= ra.rs / den;
638  }
639  }
640  }
641 
642  private:
646  const PhaseUsage phaseUsage_;
647 
651  const RegionMapping<Region> rmap_;
652 
656  struct Attributes {
657  Attributes()
658  : pressure (0.0)
659  , temperature(0.0)
660  , rs(0.0)
661  , rv(0.0)
662  , pv(0.0)
663  {}
664 
665  double pressure;
666  double temperature;
667  double rs;
668  double rv;
669  double pv;
670  };
671 
672  Details::RegionAttributes<RegionId, Attributes> attr_;
673 
674 
689  template<bool is_parallel>
690  void
691  calcAverages(const BlackoilState& state, const boost::any& info,
692  const std::vector<double>& ownerShip)
693  {
694  const auto& press = state.pressure();
695  const auto& temp = state.temperature();
696  const auto& Rv = state.rv();
697  const auto& Rs = state.gasoilratio();
698 
699  for (const auto& reg : rmap_.activeRegions()) {
700  auto& ra = attr_.attributes(reg);
701  auto& p = ra.pressure;
702  auto& T = ra.temperature;
703  auto& rs = ra.rs;
704  auto& rv = ra.rv;
705 
706  std::size_t n = 0;
707  p = T = 0.0;
708  for (const auto& cell : rmap_.cells(reg)) {
709  auto increment = Details::
710  AverageIncrementCalculator<is_parallel>()(press, temp, Rs, Rv,
711  ownerShip,
712  cell);
713  p += std::get<0>(increment);
714  T += std::get<1>(increment);
715  rs += std::get<2>(increment);
716  rv += std::get<3>(increment);
717  n += std::get<4>(increment);
718  }
719  std::size_t global_n = n;
720  double global_p = p;
721  double global_T = T;
722  double global_rs = rs;
723  double global_rv = rv;
724 #if HAVE_MPI
725  if ( is_parallel )
726  {
727  const auto& real_info = boost::any_cast<const ParallelISTLInformation&>(info);
728  global_n = real_info.communicator().sum(n);
729  global_p = real_info.communicator().sum(p);
730  global_rs = real_info.communicator().sum(rs);
731  global_rv = real_info.communicator().sum(rv);
732  global_T = real_info.communicator().sum(T);
733  }
734 #endif
735  p = global_p / global_n;
736  rs = global_rs / global_n;
737  rv = global_rv / global_n;
738  T = global_T / global_n;
739  }
740  }
741 
742 
743 
744  };
745  } // namespace RateConverter
746 } // namespace Opm
747 
748 #endif /* OPM_RATECONVERTER_HPP_HEADER_INCLUDED */
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition: RateConverter.hpp:318
Convert component rates at surface conditions to phase (voidage) rates at reservoir conditions...
Definition: RateConverter.hpp:407
const Attributes & attributes(const RegionID reg) const
Request read-only access to region&#39;s attributes.
Definition: RateConverter.hpp:209
int oil(const PhaseUsage &pu)
Numerical ID of active oil phase.
Definition: RateConverter.hpp:358
Provide mapping from Region IDs to user-specified collection of per-region attributes.
Definition: RateConverter.hpp:142
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition: RateConverter.hpp:305
std::tuple< double, double, double, double, int > operator()(const std::vector< double > &pressure, const std::vector< double > &temperature, const std::vector< double > &rs, const std::vector< double > &rv, const std::vector< double > &ownership, std::size_t cell)
Computes the temperature, pressure, and counter increment.
Definition: RateConverter.hpp:92
int cell(const RegionID reg) const
Retrieve representative cell in region.
Definition: RateConverter.hpp:196
Attributes & attributes(const RegionID reg)
Request modifiable access to region&#39;s attributes.
Definition: RateConverter.hpp:222
int water(const PhaseUsage &pu)
Numerical ID of active water phase.
Definition: RateConverter.hpp:338
RegionMapping< Region >::RegionId RegionId
Region identifier.
Definition: RateConverter.hpp:564
void defineState(const EbosSimulator &simulator)
Compute pore volume averaged hydrocarbon state pressure, rs and rv.
Definition: RateConverter.hpp:468
bool water(const PhaseUsage &pu)
Active water predicate.
Definition: RateConverter.hpp:292
int gas(const PhaseUsage &pu)
Numerical ID of active gas phase.
Definition: RateConverter.hpp:378
Computes the temperature, pressure, and counter increment.
Definition: RateConverter.hpp:78
SurfaceToReservoirVoidage(const PhaseUsage &phaseUsage, const Region &region)
Constructor.
Definition: RateConverter.hpp:416
RegionAttributes(const RMap &rmap, const Attributes &attr)
Constructor.
Definition: RateConverter.hpp:167
void calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion.
Definition: RateConverter.hpp:589
typename Select::RegionIDParameter< RegionId, std::is_integral< RegionId >::value >::type RegionID
Expose RegionId as a vocabulary type for use in query methods.
Definition: RateConverter.hpp:151
void defineState(const BlackoilState &state, const boost::any &info=boost::any())
Compute average hydrocarbon pressure and maximum dissolution and evaporation at average hydrocarbon p...
Definition: RateConverter.hpp:440