OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

split RateConverter

Browse Source
This commit is contained in:
Tor Harald Sandve
2021-09-06 14:20:40 +02:00
parent 9d2f26f7e8
commit f9832d8830
6 changed files with 704 additions and 606 deletions
Download Patch File Download Diff File Expand all files Collapse all files

636
opm/simulators/wells/RateConverter.hpp
View File

@@ -25,7 +25,7 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/grid/utility/RegionMapping.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
#include <dune/grid/common/gridenums.hh>
#include <algorithm>
#include <cmath>
@@ -47,371 +47,6 @@
namespace Opm {
namespace RateConverter {
/**
* Convenience tools for implementing the rate conversion
* facility.
*/
namespace Details {
namespace Select {
template <class RegionID, bool>
struct RegionIDParameter
{
using type =
typename std::remove_reference<RegionID>::type &;
};
template <class RegionID>
struct RegionIDParameter<RegionID, true>
{
using type = RegionID;
};
} // Select
/**
* \brief Computes the temperature, pressure, and counter increment.
*
* In a parallel run only cells owned contribute to the cell average.
* \tparam is_parallel Whether this is a parallel run.
*/
template<bool is_parallel>
struct AverageIncrementCalculator
{
/**
* \brief Computes the temperature, pressure, and counter increment.
* \param pressure The pressure.
* \param temperature The temperature.
* \param rs The rs.
* \param rv The rv.
* \param cell The current cell index.
* \param ownership A vector indicating whether a cell is owned
* by this process (value 1), or not (value 0).
* \param cell The cell index.
*/
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){
if ( ownership[cell] )
{
return std::make_tuple(pressure[cell],
temperature[cell],
rs[cell],
rv[cell],
1);
}
else
{
return std::make_tuple(0, 0, 0, 0, 0);
}
}
};
template<>
struct AverageIncrementCalculator<false>
{
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>&,
std::size_t cell){
return std::make_tuple(pressure[cell],
temperature[cell],
rs[cell],
rv[cell],
1);
}
};
/**
* Provide mapping from Region IDs to user-specified collection
* of per-region attributes.
*
* \tparam RegionId Region identifier type. Must be hashable by
* \code std::hash<> \endcode. Typically a built-in integer
* type--e.g., \c int.
*
* \tparam Attributes User-defined type that represents
* collection of attributes that have meaning in a per-region
* aggregate sense. Must be copy-constructible.
*/
template <typename RegionId, class Attributes>
class RegionAttributes
{
public:
/**
* Expose \c RegionId as a vocabulary type for use in query
* methods.
*/
using RegionID =
typename Select::RegionIDParameter
<RegionId, std::is_integral<RegionId>::value>::type;
/**
* Constructor.
*
* \tparam RMap Class type that implements the RegionMapping
* protocol. Typically an instantiation of \code
* Opm::RegionMapping<> \endcode.
*
* \param[in] rmap Specific region mapping that provides
* reverse lookup from regions to cells.
*
* \param[in] attr Pre-constructed initialiser for \c
* Attributes.
*/
template <class RMap>
RegionAttributes(const RMap& rmap,
const Attributes& attr)
{
using VT = typename AttributeMap::value_type;
for (const auto& r : rmap.activeRegions()) {
auto v = std::make_unique<Value>(attr);
const auto stat = attr_.insert(VT(r, std::move(v)));
if (stat.second) {
// New value inserted.
const auto& cells = rmap.cells(r);
assert (! cells.empty());
// Region's representative cell.
stat.first->second->cell_ = cells[0];
}
}
}
/**
* Retrieve representative cell in region.
*
* \param[in] reg Specific region.
*
* \return Representative cell in region \p reg.
*/
int cell(const RegionID reg) const
{
return this->find(reg).cell_;
}
bool has(const RegionID reg) const
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
return false;
}
return true;
}
void insert(const RegionID r, const Attributes& attr)
{
using VT = typename AttributeMap::value_type;
auto v = std::make_unique<Value>(attr);
const auto stat = attr_.insert(VT(r, std::move(v)));
if (stat.second) {
// Region's representative cell.
stat.first->second->cell_ = -1.0;
}
}
/**
* Request read-only access to region's attributes.
*
* \param[in] reg Specific region.
*
* \return Read-only access to region \p reg's per-region
* attributes.
*/
const Attributes& attributes(const RegionID reg) const
{
return this->find(reg).attr_;
}
/**
* Request modifiable access to region's attributes.
*
* \param[in] reg Specific region.
*
* \return Read-write access to region \p reg's per-region
* attributes.
*/
Attributes& attributes(const RegionID reg)
{
return this->find(reg).attr_;
}
private:
/**
* Aggregate per-region attributes along with region's
* representative cell.
*/
struct Value {
Value(const Attributes& attr)
: attr_(attr)
, cell_(-1)
{}
Attributes attr_;
int cell_;
};
using ID =
typename std::remove_reference<RegionId>::type;
using AttributeMap =
std::unordered_map<ID, std::unique_ptr<Value>>;
AttributeMap attr_;
/**
* Read-only access to region's properties.
*/
const Value& find(const RegionID reg) const
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
throw std::invalid_argument("Unknown region ID");
}
return *i->second;
}
/**
* Read-write access to region's properties.
*/
Value& find(const RegionID reg)
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
throw std::invalid_argument("Unknown region ID");
}
return *i->second;
}
};
/**
* Convenience functions for querying presence/absence of
* active phases.
*/
namespace PhaseUsed {
/**
* Active water predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not water is an active phase.
*/
inline bool
water(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
}
/**
* Active oil predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not oil is an active phase.
*/
inline bool
oil(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
}
/**
* Active gas predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not gas is an active phase.
*/
inline bool
gas(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
}
} // namespace PhaseUsed
/**
* Convenience functions for querying numerical IDs
* ("positions") of active phases.
*/
namespace PhasePos {
/**
* Numerical ID of active water phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of water if
* active, -1 if not.
*/
inline int
water(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::water(pu)) {
p = pu.phase_pos[ BlackoilPhases::Aqua ];
}
return p;
}
/**
* Numerical ID of active oil phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of oil if
* active, -1 if not.
*/
inline int
oil(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::oil(pu)) {
p = pu.phase_pos[ BlackoilPhases::Liquid ];
}
return p;
}
/**
* Numerical ID of active gas phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of gas if
* active, -1 if not.
*/
inline int
gas(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::gas(pu)) {
p = pu.phase_pos[ BlackoilPhases::Vapour ];
}
return p;
}
} // namespace PhasePos
} // namespace Details
/**
* Convert component rates at surface conditions to phase
@@ -510,7 +145,7 @@ namespace Opm {
// only count oil and gas filled parts of the domain
double hydrocarbon = 1.0;
const auto& pu = phaseUsage_;
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
}
@@ -522,15 +157,15 @@ namespace Opm {
if (hydrocarbonPV > 0.) {
auto& attr = attributes_hpv[reg];
attr.pv += hydrocarbonPV;
if (Details::PhaseUsed::oil(pu) && Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu) && RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.rs += fs.Rs().value() * hydrocarbonPV;
attr.rv += fs.Rv().value() * hydrocarbonPV;
}
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
} else {
assert(Details::PhaseUsed::gas(pu));
assert(RegionAttributeHelpers::PhaseUsed::gas(pu));
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
}
@@ -540,18 +175,18 @@ namespace Opm {
if (pv_cell > 0.) {
auto& attr = attributes_pv[reg];
attr.pv += pv_cell;
if (Details::PhaseUsed::oil(pu) && Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu) && RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.rs += fs.Rs().value() * pv_cell;
attr.rv += fs.Rv().value() * pv_cell;
}
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell;
} else if (Details::PhaseUsed::gas(pu)) {
} else if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * pv_cell;
} else {
assert(Details::PhaseUsed::water(pu));
assert(RegionAttributeHelpers::PhaseUsed::water(pu));
attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell;
}
@@ -642,13 +277,13 @@ namespace Opm {
const double T = ra.temperature;
const double saltConcentration = ra.saltConcentration;
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
const int iw = RegionAttributeHelpers::PhasePos::water(pu);
const int io = RegionAttributeHelpers::PhasePos::oil (pu);
const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
@@ -663,7 +298,7 @@ namespace Opm {
// Determinant of 'R' matrix
const double detR = 1.0 - (Rs * Rv);
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
@@ -671,12 +306,12 @@ namespace Opm {
coeff[io] += 1.0 / den;
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
coeff[ig] -= ra.rv / den;
}
}
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
@@ -684,7 +319,7 @@ namespace Opm {
coeff[ig] += 1.0 / den;
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
coeff[io] -= ra.rs / den;
}
}
@@ -700,13 +335,13 @@ namespace Opm {
const double T = ra.temperature;
const double saltConcentration = ra.saltConcentration;
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
const int iw = RegionAttributeHelpers::PhasePos::water(pu);
const int io = RegionAttributeHelpers::PhasePos::oil (pu);
const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
@@ -714,12 +349,12 @@ namespace Opm {
coeff[iw] = 1.0 / bw;
}
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
coeff[io] += 1.0 / bo;
}
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
coeff[ig] += 1.0 / bg;
}
@@ -758,11 +393,11 @@ namespace Opm {
const double T = ra.temperature;
const double saltConcentration = ra.saltConcentration;
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
const int iw = RegionAttributeHelpers::PhasePos::water(pu);
const int io = RegionAttributeHelpers::PhasePos::oil (pu);
const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
@@ -790,7 +425,7 @@ namespace Opm {
// Determinant of 'R' matrix
const double detR = 1.0 - (Rs * Rv);
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
@@ -798,14 +433,14 @@ namespace Opm {
voidage_rates[io] = surface_rates[io];
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
voidage_rates[io] -= Rv * surface_rates[ig];
}
voidage_rates[io] /= den;
}
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
@@ -813,7 +448,7 @@ namespace Opm {
voidage_rates[ig] = surface_rates[ig];
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
voidage_rates[ig] -= Rs * surface_rates[io];
}
@@ -878,214 +513,7 @@ namespace Opm {
double saltConcentration;
};
Details::RegionAttributes<RegionId, Attributes> attr_;
};
/**
* Computes hydrocarbon weighed average pressures over regions
*
* \tparam FluidSystem Fluid system class. Expected to be a BlackOilFluidSystem
*
* \tparam Region Type of a forward region mapping. Expected
* to provide indexed access through \code operator[]()
* \endcode as well as inner types \c value_type, \c
* size_type, and \c const_iterator. Typically \code
* std::vector<int> \endcode.
*/
template <class FluidSystem, class Region>
class AverageRegionalPressure {
public:
/**
* Constructor.
*
* \param[in] region Forward region mapping. Often
* corresponds to the "FIPNUM" mapping of an ECLIPSE input
* deck.
*/
AverageRegionalPressure(const PhaseUsage& phaseUsage,
const Region& region)
: phaseUsage_(phaseUsage)
, rmap_ (region)
, attr_ (rmap_, Attributes())
{
}
/**
* Compute pore volume averaged hydrocarbon state pressure, *
*/
template <typename ElementContext, class EbosSimulator>
void defineState(const EbosSimulator& simulator)
{
int numRegions = 0;
const auto& gridView = simulator.gridView();
const auto& comm = gridView.comm();
for (const auto& reg : rmap_.activeRegions()) {
numRegions = std::max(numRegions, reg);
}
numRegions = comm.max(numRegions);
for (int reg = 0; reg < numRegions; ++ reg) {
if(!attr_.has(reg))
attr_.insert(reg, Attributes());
}
// create map from cell to region
// and set all attributes to zero
for (int reg = 0; reg < numRegions; ++ reg) {
auto& ra = attr_.attributes(reg);
ra.pressure = 0.0;
ra.pv = 0.0;
}
// quantities for pore volume average
std::unordered_map<RegionId, Attributes> attributes_pv;
// quantities for hydrocarbon volume average
std::unordered_map<RegionId, Attributes> attributes_hpv;
for (int reg = 0; reg < numRegions; ++ reg) {
attributes_pv.insert({reg, Attributes()});
attributes_hpv.insert({reg, Attributes()});
}
for (int reg = 0; reg < numRegions; ++ reg) {
auto& ra = attributes_pv[reg];
ra.pressure = 0.0;
ra.pv = 0.0;
}
for (int reg = 0; reg < numRegions; ++ reg) {
auto& ra = attributes_hpv[reg];
ra.pressure = 0.0;
ra.pv = 0.0;
}
ElementContext elemCtx( simulator );
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// use pore volume weighted averages.
const double pv_cell =
simulator.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
// only count oil and gas filled parts of the domain
double hydrocarbon = 1.0;
const auto& pu = phaseUsage_;
if (Details::PhaseUsed::water(pu)) {
hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
}
const int reg = rmap_.region(cellIdx);
assert(reg >= 0);
// sum p, rs, rv, and T.
const double hydrocarbonPV = pv_cell*hydrocarbon;
if (hydrocarbonPV > 0.) {
auto& attr = attributes_hpv[reg];
attr.pv += hydrocarbonPV;
if (Details::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
} else {
assert(Details::PhaseUsed::gas(pu));
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
}
}
if (pv_cell > 0.) {
auto& attr = attributes_pv[reg];
attr.pv += pv_cell;
if (Details::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
} else if (Details::PhaseUsed::gas(pu)) {
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
} else {
assert(Details::PhaseUsed::water(pu));
attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
}
}
}
for (int reg = 0; reg < numRegions; ++ reg) {
auto& ra = attr_.attributes(reg);
const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
// TODO: should we have some epsilon here instead of zero?
if (hpv_sum > 0.) {
const auto& attri_hpv = attributes_hpv[reg];
const double p_hpv_sum = comm.sum(attri_hpv.pressure);
ra.pressure = p_hpv_sum / hpv_sum;
} else {
// using the pore volume to do the averaging
const auto& attri_pv = attributes_pv[reg];
const double pv_sum = comm.sum(attri_pv.pv);
assert(pv_sum > 0.);
const double p_pv_sum = comm.sum(attri_pv.pressure);
ra.pressure = p_pv_sum / pv_sum;
}
}
}
/**
* Region identifier.
*
* Integral type.
*/
typedef typename RegionMapping<Region>::RegionId RegionId;
/**
* Average pressure
*
*/
double
fpr(const RegionId r) const
{
const auto& ra = attr_.attributes(r);
return ra.pressure;
}
private:
/**
* Fluid property object.
*/
const PhaseUsage phaseUsage_;
/**
* "Fluid-in-place" region mapping (forward and reverse).
*/
const RegionMapping<Region> rmap_;
/**
* Derived property attributes for each active region.
*/
struct Attributes {
Attributes()
: pressure (0.0)
, pv(0.0)
{}
double pressure;
double pv;
};
Details::RegionAttributes<RegionId, Attributes> attr_;
RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
};
} // namespace RateConverter