21 #ifndef OPM_RATECONVERTER_HPP_HEADER_INCLUDED 22 #define OPM_RATECONVERTER_HPP_HEADER_INCLUDED 24 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp> 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> 35 #include <type_traits> 36 #include <unordered_map> 50 namespace RateConverter {
57 template <
class RegionID,
bool>
61 typename std::remove_reference<RegionID>::type &;
64 template <
class RegionID>
67 using type = RegionID;
77 template<
bool is_parallel>
91 std::tuple<double, double, double, double, int>
93 const std::vector<double>& temperature,
94 const std::vector<double>& rs,
95 const std::vector<double>& rv,
96 const std::vector<double>& ownership,
98 if ( ownership[cell] )
100 return std::make_tuple(pressure[cell],
108 return std::make_tuple(0, 0, 0, 0, 0);
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>&,
122 return std::make_tuple(pressure[cell],
141 template <
typename RegionId,
class Attributes>
151 <RegionId, std::is_integral<RegionId>::value>::type;
166 template <
class RMap>
168 const Attributes& attr)
170 using VT =
typename AttributeMap::value_type;
172 for (
const auto& r : rmap.activeRegions()) {
173 auto v = std::unique_ptr<Value>(
new Value(attr));
175 const auto stat = attr_.insert(VT(r, std::move(v)));
179 const auto& cells = rmap.cells(r);
181 assert (! cells.empty());
184 stat.first->second->cell_ = cells[0];
198 return this->find(reg).cell_;
211 return this->find(reg).attr_;
224 return this->find(reg).attr_;
233 Value(
const Attributes& attr)
243 typename std::remove_reference<RegionId>::type;
246 std::unordered_map<ID, std::unique_ptr<Value>>;
253 const Value& find(
const RegionID reg)
const 255 const auto& i = attr_.find(reg);
257 if (i == attr_.end()) {
258 throw std::invalid_argument(
"Unknown region ID");
269 const auto& i = attr_.find(reg);
271 if (i == attr_.end()) {
272 throw std::invalid_argument(
"Unknown region ID");
283 namespace PhaseUsed {
294 return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
307 return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
320 return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
343 p = pu.phase_pos[ BlackoilPhases::Aqua ];
363 p = pu.phase_pos[ BlackoilPhases::Liquid ];
383 p = pu.phase_pos[ BlackoilPhases::Vapour ];
406 template <
class Flu
idSystem,
class Region>
417 const Region& region)
418 : phaseUsage_(phaseUsage)
420 , attr_ (rmap_, Attributes())
441 const boost::any& info = boost::any())
444 if( info.type() ==
typeid(ParallelISTLInformation) )
446 const auto& ownership =
447 boost::any_cast<
const ParallelISTLInformation&>(info)
448 .updateOwnerMask(state.pressure());
449 calcAverages<true>(state, info, ownership);
454 std::vector<double> dummyOwnership;
455 calcAverages<false>(state, info, dummyOwnership);
467 template <
typename ElementContext,
class EbosSimulator>
473 const auto& grid = simulator.gridManager().grid();
474 const unsigned numCells = grid.size(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;
480 auto& ra = attr_.attributes(reg);
482 ra.temperature = 0.0;
489 ElementContext elemCtx( simulator );
490 const auto& gridView = simulator.gridView();
491 const auto& comm = gridView.comm();
493 const auto& elemEndIt = gridView.template end<0>();
494 for (
auto elemIt = gridView.template begin</*codim=*/0>();
499 const auto& elem = *elemIt;
500 if (elem.partitionType() != Dune::InteriorEntity)
503 elemCtx.updatePrimaryStencil(elem);
504 elemCtx.updatePrimaryIntensiveQuantities(0);
505 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
506 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
507 const auto& fs = intQuants.fluidState();
509 const double pv_cell =
510 simulator.model().dofTotalVolume(cellIdx)
511 * intQuants.porosity().value();
514 double hydrocarbon = 1.0;
515 const auto& pu = phaseUsage_;
517 hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
520 int reg = cell2region[cellIdx];
522 auto& ra = attr_.attributes(reg);
523 auto& p = ra.pressure;
524 auto& T = ra.temperature;
530 double hydrocarbonPV = pv_cell*hydrocarbon;
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;
538 for (
const auto& reg : rmap_.activeRegions()) {
539 auto& ra = attr_.attributes(reg);
540 auto& p = ra.pressure;
541 auto& T = ra.temperature;
564 typedef typename RegionMapping<Region>::RegionId
RegionId;
587 template <
class Coeff>
591 const auto& pu = phaseUsage_;
592 const auto& ra = attr_.attributes(r);
593 const double p = ra.pressure;
594 const double T = ra.temperature;
600 std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
605 const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p);
607 coeff[iw] = 1.0 / bw;
611 const double detR = 1.0 - (ra.rs * ra.rv);
616 const double Rs = ra.rs;
617 const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
618 const double den = bo * detR;
620 coeff[io] += 1.0 / den;
623 coeff[ig] -= ra.rv / den;
630 const double Rv = ra.rv;
631 const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
632 const double den = bg * detR;
634 coeff[ig] += 1.0 / den;
637 coeff[io] -= ra.rs / den;
646 const PhaseUsage phaseUsage_;
651 const RegionMapping<Region> rmap_;
672 Details::RegionAttributes<RegionId, Attributes> attr_;
689 template<
bool is_parallel>
691 calcAverages(
const BlackoilState& state,
const boost::any& info,
692 const std::vector<double>& ownerShip)
694 const auto& press = state.pressure();
695 const auto& temp = state.temperature();
696 const auto& Rv = state.rv();
697 const auto& Rs = state.gasoilratio();
699 for (
const auto& reg : rmap_.activeRegions()) {
700 auto& ra = attr_.attributes(reg);
701 auto& p = ra.pressure;
702 auto& T = ra.temperature;
708 for (
const auto& cell : rmap_.cells(reg)) {
709 auto increment = Details::
710 AverageIncrementCalculator<is_parallel>()(press, temp, Rs, Rv,
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);
719 std::size_t global_n = n;
722 double global_rs = rs;
723 double global_rv = rv;
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);
735 p = global_p / global_n;
736 rs = global_rs / global_n;
737 rv = global_rv / global_n;
738 T = global_T / global_n;
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's attributes.
Definition: RateConverter.hpp:209
int oil(const PhaseUsage &pu)
Numerical ID of active oil phase.
Definition: RateConverter.hpp:358
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
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
void calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion.
Definition: RateConverter.hpp:589
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
Attributes & attributes(const RegionID reg)
Request modifiable access to region'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 ®ion)
Constructor.
Definition: RateConverter.hpp:416
RegionAttributes(const RMap &rmap, const Attributes &attr)
Constructor.
Definition: RateConverter.hpp:167
Definition: RateConverter.hpp:58
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
int cell(const RegionID reg) const
Retrieve representative cell in region.
Definition: RateConverter.hpp:196
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