00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED
00022 #define OPM_RATECONVERTER_HPP_HEADER_INCLUDED
00023
00024 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
00025
00026 #include <opm/core/props/BlackoilPhases.hpp>
00027 #include <opm/core/simulator/BlackoilState.hpp>
00028 #include <opm/core/utility/RegionMapping.hpp>
00029 #include <opm/core/linalg/ParallelIstlInformation.hpp>
00030
00031 #include <algorithm>
00032 #include <cmath>
00033 #include <memory>
00034 #include <stdexcept>
00035 #include <type_traits>
00036 #include <unordered_map>
00037 #include <utility>
00038 #include <vector>
00049 namespace Opm {
00050 namespace RateConverter {
00055 namespace Details {
00056 namespace Select {
00057 template <class RegionID, bool>
00058 struct RegionIDParameter
00059 {
00060 using type =
00061 typename std::remove_reference<RegionID>::type &;
00062 };
00063
00064 template <class RegionID>
00065 struct RegionIDParameter<RegionID, true>
00066 {
00067 using type = RegionID;
00068 };
00069 }
00070
00077 template<bool is_parallel>
00078 struct AverageIncrementCalculator
00079 {
00091 std::tuple<double, double, double, double, int>
00092 operator()(const std::vector<double>& pressure,
00093 const std::vector<double>& temperature,
00094 const std::vector<double>& rs,
00095 const std::vector<double>& rv,
00096 const std::vector<double>& ownership,
00097 std::size_t cell){
00098 if ( ownership[cell] )
00099 {
00100 return std::make_tuple(pressure[cell],
00101 temperature[cell],
00102 rs[cell],
00103 rv[cell],
00104 1);
00105 }
00106 else
00107 {
00108 return std::make_tuple(0, 0, 0, 0, 0);
00109 }
00110 }
00111 };
00112 template<>
00113 struct AverageIncrementCalculator<false>
00114 {
00115 std::tuple<double, double, double, double, int>
00116 operator()(const std::vector<double>& pressure,
00117 const std::vector<double>& temperature,
00118 const std::vector<double>& rs,
00119 const std::vector<double>& rv,
00120 const std::vector<double>&,
00121 std::size_t cell){
00122 return std::make_tuple(pressure[cell],
00123 temperature[cell],
00124 rs[cell],
00125 rv[cell],
00126 1);
00127 }
00128 };
00141 template <typename RegionId, class Attributes>
00142 class RegionAttributes
00143 {
00144 public:
00149 using RegionID =
00150 typename Select::RegionIDParameter
00151 <RegionId, std::is_integral<RegionId>::value>::type;
00152
00166 template <class RMap>
00167 RegionAttributes(const RMap& rmap,
00168 const Attributes& attr)
00169 {
00170 using VT = typename AttributeMap::value_type;
00171
00172 for (const auto& r : rmap.activeRegions()) {
00173 auto v = std::unique_ptr<Value>(new Value(attr));
00174
00175 const auto stat = attr_.insert(VT(r, std::move(v)));
00176
00177 if (stat.second) {
00178
00179 const auto& cells = rmap.cells(r);
00180
00181 assert (! cells.empty());
00182
00183
00184 stat.first->second->cell_ = cells[0];
00185 }
00186 }
00187 }
00188
00196 int cell(const RegionID reg) const
00197 {
00198 return this->find(reg).cell_;
00199 }
00200
00209 const Attributes& attributes(const RegionID reg) const
00210 {
00211 return this->find(reg).attr_;
00212 }
00213
00222 Attributes& attributes(const RegionID reg)
00223 {
00224 return this->find(reg).attr_;
00225 }
00226
00227 private:
00232 struct Value {
00233 Value(const Attributes& attr)
00234 : attr_(attr)
00235 , cell_(-1)
00236 {}
00237
00238 Attributes attr_;
00239 int cell_;
00240 };
00241
00242 using ID =
00243 typename std::remove_reference<RegionId>::type;
00244
00245 using AttributeMap =
00246 std::unordered_map<ID, std::unique_ptr<Value>>;
00247
00248 AttributeMap attr_;
00249
00253 const Value& find(const RegionID reg) const
00254 {
00255 const auto& i = attr_.find(reg);
00256
00257 if (i == attr_.end()) {
00258 throw std::invalid_argument("Unknown region ID");
00259 }
00260
00261 return *i->second;
00262 }
00263
00267 Value& find(const RegionID reg)
00268 {
00269 const auto& i = attr_.find(reg);
00270
00271 if (i == attr_.end()) {
00272 throw std::invalid_argument("Unknown region ID");
00273 }
00274
00275 return *i->second;
00276 }
00277 };
00278
00283 namespace PhaseUsed {
00291 inline bool
00292 water(const PhaseUsage& pu)
00293 {
00294 return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
00295 }
00296
00304 inline bool
00305 oil(const PhaseUsage& pu)
00306 {
00307 return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
00308 }
00309
00317 inline bool
00318 gas(const PhaseUsage& pu)
00319 {
00320 return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
00321 }
00322 }
00323
00328 namespace PhasePos {
00337 inline int
00338 water(const PhaseUsage& pu)
00339 {
00340 int p = -1;
00341
00342 if (PhaseUsed::water(pu)) {
00343 p = pu.phase_pos[ BlackoilPhases::Aqua ];
00344 }
00345
00346 return p;
00347 }
00348
00357 inline int
00358 oil(const PhaseUsage& pu)
00359 {
00360 int p = -1;
00361
00362 if (PhaseUsed::oil(pu)) {
00363 p = pu.phase_pos[ BlackoilPhases::Liquid ];
00364 }
00365
00366 return p;
00367 }
00368
00377 inline int
00378 gas(const PhaseUsage& pu)
00379 {
00380 int p = -1;
00381
00382 if (PhaseUsed::gas(pu)) {
00383 p = pu.phase_pos[ BlackoilPhases::Vapour ];
00384 }
00385
00386 return p;
00387 }
00388 }
00389 }
00390
00406 template <class FluidSystem, class Region>
00407 class SurfaceToReservoirVoidage {
00408 public:
00416 SurfaceToReservoirVoidage(const PhaseUsage& phaseUsage,
00417 const Region& region)
00418 : phaseUsage_(phaseUsage)
00419 , rmap_ (region)
00420 , attr_ (rmap_, Attributes())
00421 {
00422 }
00423
00439 void
00440 defineState(const BlackoilState& state,
00441 const boost::any& info = boost::any())
00442 {
00443 #if HAVE_MPI
00444 if( info.type() == typeid(ParallelISTLInformation) )
00445 {
00446 const auto& ownership =
00447 boost::any_cast<const ParallelISTLInformation&>(info)
00448 .updateOwnerMask(state.pressure());
00449 calcAverages<true>(state, info, ownership);
00450 }
00451 else
00452 #endif
00453 {
00454 std::vector<double> dummyOwnership;
00455 calcAverages<false>(state, info, dummyOwnership);
00456 }
00457 }
00458
00467 template <typename ElementContext, class EbosSimulator>
00468 void defineState(const EbosSimulator& simulator)
00469 {
00470
00471
00472
00473 const auto& grid = simulator.gridManager().grid();
00474 const unsigned numCells = grid.size(0);
00475 std::vector<int> cell2region(numCells, -1);
00476 for (const auto& reg : rmap_.activeRegions()) {
00477 for (const auto& cell : rmap_.cells(reg)) {
00478 cell2region[cell] = reg;
00479 }
00480 auto& ra = attr_.attributes(reg);
00481 ra.pressure = 0.0;
00482 ra.temperature = 0.0;
00483 ra.rs = 0.0;
00484 ra.rv = 0.0;
00485 ra.pv = 0.0;
00486
00487 }
00488
00489 ElementContext elemCtx( simulator );
00490 const auto& gridView = simulator.gridView();
00491 const auto& comm = gridView.comm();
00492
00493 const auto& elemEndIt = gridView.template end<0>();
00494 for (auto elemIt = gridView.template begin</*codim=*/0>();
00495 elemIt != elemEndIt;
00496 ++elemIt)
00497 {
00498
00499 const auto& elem = *elemIt;
00500 if (elem.partitionType() != Dune::InteriorEntity)
00501 continue;
00502
00503 elemCtx.updatePrimaryStencil(elem);
00504 elemCtx.updatePrimaryIntensiveQuantities(0);
00505 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
00506 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
00507 const auto& fs = intQuants.fluidState();
00508
00509 const double pv_cell =
00510 simulator.model().dofTotalVolume(cellIdx)
00511 * intQuants.porosity().value();
00512
00513
00514 double hydrocarbon = 1.0;
00515 const auto& pu = phaseUsage_;
00516 if (Details::PhaseUsed::water(pu)) {
00517 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
00518 }
00519
00520 int reg = cell2region[cellIdx];
00521 assert(reg >= 0);
00522 auto& ra = attr_.attributes(reg);
00523 auto& p = ra.pressure;
00524 auto& T = ra.temperature;
00525 auto& rs = ra.rs;
00526 auto& rv = ra.rv;
00527 auto& pv = ra.pv;
00528
00529
00530 double hydrocarbonPV = pv_cell*hydrocarbon;
00531 pv += hydrocarbonPV;
00532 p += fs.pressure(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
00533 rs += fs.Rs().value()*hydrocarbonPV;
00534 rv += fs.Rv().value()*hydrocarbonPV;
00535 T += fs.temperature(FluidSystem::oilPhaseIdx).value()*hydrocarbonPV;
00536 }
00537
00538 for (const auto& reg : rmap_.activeRegions()) {
00539 auto& ra = attr_.attributes(reg);
00540 auto& p = ra.pressure;
00541 auto& T = ra.temperature;
00542 auto& rs = ra.rs;
00543 auto& rv = ra.rv;
00544 auto& pv = ra.pv;
00545
00546 p = comm.sum(p);
00547 T = comm.sum(T);
00548 rs = comm.sum(rs);
00549 rv = comm.sum(rv);
00550 pv = comm.sum(pv);
00551
00552 p /= pv;
00553 T /= pv;
00554 rs /= pv;
00555 rv /= pv;
00556 }
00557 }
00558
00564 typedef typename RegionMapping<Region>::RegionId RegionId;
00565
00587 template <class Coeff>
00588 void
00589 calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const
00590 {
00591 const auto& pu = phaseUsage_;
00592 const auto& ra = attr_.attributes(r);
00593 const double p = ra.pressure;
00594 const double T = ra.temperature;
00595
00596 const int iw = Details::PhasePos::water(pu);
00597 const int io = Details::PhasePos::oil (pu);
00598 const int ig = Details::PhasePos::gas (pu);
00599
00600 std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
00601
00602 if (Details::PhaseUsed::water(pu)) {
00603
00604
00605 const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
00606
00607 coeff[iw] = 1.0 / bw;
00608 }
00609
00610
00611 const double detR = 1.0 - (ra.rs * ra.rv);
00612
00613 if (Details::PhaseUsed::oil(pu)) {
00614
00615
00616 const double Rs = ra.rs;
00617 const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
00618 const double den = bo * detR;
00619
00620 coeff[io] += 1.0 / den;
00621
00622 if (Details::PhaseUsed::gas(pu)) {
00623 coeff[ig] -= ra.rv / den;
00624 }
00625 }
00626
00627 if (Details::PhaseUsed::gas(pu)) {
00628
00629
00630 const double Rv = ra.rv;
00631 const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
00632 const double den = bg * detR;
00633
00634 coeff[ig] += 1.0 / den;
00635
00636 if (Details::PhaseUsed::oil(pu)) {
00637 coeff[io] -= ra.rs / den;
00638 }
00639 }
00640 }
00641
00642 private:
00646 const PhaseUsage phaseUsage_;
00647
00651 const RegionMapping<Region> rmap_;
00652
00656 struct Attributes {
00657 Attributes()
00658 : pressure (0.0)
00659 , temperature(0.0)
00660 , rs(0.0)
00661 , rv(0.0)
00662 , pv(0.0)
00663 {}
00664
00665 double pressure;
00666 double temperature;
00667 double rs;
00668 double rv;
00669 double pv;
00670 };
00671
00672 Details::RegionAttributes<RegionId, Attributes> attr_;
00673
00674
00689 template<bool is_parallel>
00690 void
00691 calcAverages(const BlackoilState& state, const boost::any& info,
00692 const std::vector<double>& ownerShip)
00693 {
00694 const auto& press = state.pressure();
00695 const auto& temp = state.temperature();
00696 const auto& Rv = state.rv();
00697 const auto& Rs = state.gasoilratio();
00698
00699 for (const auto& reg : rmap_.activeRegions()) {
00700 auto& ra = attr_.attributes(reg);
00701 auto& p = ra.pressure;
00702 auto& T = ra.temperature;
00703 auto& rs = ra.rs;
00704 auto& rv = ra.rv;
00705
00706 std::size_t n = 0;
00707 p = T = 0.0;
00708 for (const auto& cell : rmap_.cells(reg)) {
00709 auto increment = Details::
00710 AverageIncrementCalculator<is_parallel>()(press, temp, Rs, Rv,
00711 ownerShip,
00712 cell);
00713 p += std::get<0>(increment);
00714 T += std::get<1>(increment);
00715 rs += std::get<2>(increment);
00716 rv += std::get<3>(increment);
00717 n += std::get<4>(increment);
00718 }
00719 std::size_t global_n = n;
00720 double global_p = p;
00721 double global_T = T;
00722 double global_rs = rs;
00723 double global_rv = rv;
00724 #if HAVE_MPI
00725 if ( is_parallel )
00726 {
00727 const auto& real_info = boost::any_cast<const ParallelISTLInformation&>(info);
00728 global_n = real_info.communicator().sum(n);
00729 global_p = real_info.communicator().sum(p);
00730 global_rs = real_info.communicator().sum(rs);
00731 global_rv = real_info.communicator().sum(rv);
00732 global_T = real_info.communicator().sum(T);
00733 }
00734 #endif
00735 p = global_p / global_n;
00736 rs = global_rs / global_n;
00737 rv = global_rv / global_n;
00738 T = global_T / global_n;
00739 }
00740 }
00741
00742
00743
00744 };
00745 }
00746 }
00747
00748 #endif