diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4ad30d701..b0da7681e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -309,6 +309,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/PerfData.hpp opm/simulators/wells/PerforationData.hpp opm/simulators/wells/RateConverter.hpp + opm/simulators/wells/RegionAttributeHelpers.hpp + opm/simulators/wells/RegionAverageCalculator.hpp opm/simulators/utils/readDeck.hpp opm/simulators/wells/SingleWellState.hpp opm/simulators/wells/TargetCalculator.hpp diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index e685ceebe..83bab1b1a 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -60,6 +60,7 @@ #include #include #include +#include #include #include #include @@ -133,7 +134,7 @@ namespace Opm { SurfaceToReservoirVoidage >; // For computing average pressured used by gpmaint - using AverageRegionalPressureType = RateConverter:: + using AverageRegionalPressureType = RegionAverageCalculator:: AverageRegionalPressure >; BlackoilWellModel(Simulator& ebosSimulator); diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index ab25d164d..2664d81ed 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -25,7 +25,7 @@ #include #include #include - +#include #include #include #include @@ -47,371 +47,6 @@ namespace Opm { namespace RateConverter { - /** - * Convenience tools for implementing the rate conversion - * facility. - */ - namespace Details { - namespace Select { - template - struct RegionIDParameter - { - using type = - typename std::remove_reference::type &; - }; - - template - struct RegionIDParameter - { - 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 - 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 - operator()(const std::vector& pressure, - const std::vector& temperature, - const std::vector& rs, - const std::vector& rv, - const std::vector& 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 - { - std::tuple - operator()(const std::vector& pressure, - const std::vector& temperature, - const std::vector& rs, - const std::vector& rv, - const std::vector&, - 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 - class RegionAttributes - { - public: - /** - * Expose \c RegionId as a vocabulary type for use in query - * methods. - */ - using RegionID = - typename Select::RegionIDParameter - ::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 - RegionAttributes(const RMap& rmap, - const Attributes& attr) - { - using VT = typename AttributeMap::value_type; - - for (const auto& r : rmap.activeRegions()) { - auto v = std::make_unique(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(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::type; - - using AttributeMap = - std::unordered_map>; - - 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 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 \endcode. - */ - template - 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 - 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 attributes_pv; - - // quantities for hydrocarbon volume average - std::unordered_map 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(); - for (auto elemIt = gridView.template begin(); - 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::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 rmap_; - - /** - * Derived property attributes for each active region. - */ - struct Attributes { - Attributes() - : pressure (0.0) - , pv(0.0) - - {} - - double pressure; - double pv; - - }; - - Details::RegionAttributes attr_; + RegionAttributeHelpers::RegionAttributes attr_; }; } // namespace RateConverter diff --git a/opm/simulators/wells/RegionAttributeHelpers.hpp b/opm/simulators/wells/RegionAttributeHelpers.hpp new file mode 100644 index 000000000..e65e7a32e --- /dev/null +++ b/opm/simulators/wells/RegionAttributeHelpers.hpp @@ -0,0 +1,409 @@ +/* + Copyright 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Statoil ASA. + Copyright 2017, IRIS + Copyright 2017, Equinor + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED +#define OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + namespace RegionAttributeHelpers { + /** + * Convenience tools for processing region + * spesific attributes + */ + namespace Select { + template + struct RegionIDParameter + { + using type = + typename std::remove_reference::type &; + }; + + template + struct RegionIDParameter + { + 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 + 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 + operator()(const std::vector& pressure, + const std::vector& temperature, + const std::vector& rs, + const std::vector& rv, + const std::vector& 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 + { + std::tuple + operator()(const std::vector& pressure, + const std::vector& temperature, + const std::vector& rs, + const std::vector& rv, + const std::vector&, + 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 + class RegionAttributes + { + public: + /** + * Expose \c RegionId as a vocabulary type for use in query + * methods. + */ + using RegionID = + typename Select::RegionIDParameter + ::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 + RegionAttributes(const RMap& rmap, + const Attributes& attr) + { + using VT = typename AttributeMap::value_type; + + for (const auto& r : rmap.activeRegions()) { + auto v = std::make_unique(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(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::type; + + using AttributeMap = + std::unordered_map>; + + 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 RegionAttributesHelpers +} // namespace Opm + +#endif /* OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED */ diff --git a/opm/simulators/wells/RegionAverageCalculator.hpp b/opm/simulators/wells/RegionAverageCalculator.hpp new file mode 100644 index 000000000..df7663944 --- /dev/null +++ b/opm/simulators/wells/RegionAverageCalculator.hpp @@ -0,0 +1,258 @@ +/* + Copyright 2021, Equinor + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED +#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +/** + * \file + * Facility for converting component rates at surface conditions to + * phase (voidage) rates at reservoir conditions. + * + * This uses the average hydrocarbon pressure to define fluid + * properties. The facility is intended to support Reservoir Voidage + * rates only ('RESV'). + */ + +namespace Opm { + namespace RegionAverageCalculator { + + /** + * 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 \endcode. + */ + template + 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 + 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 attributes_pv; + + // quantities for hydrocarbon volume average + std::unordered_map 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(); + for (auto elemIt = gridView.template begin(); + 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 (RegionAttributeHelpers::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 (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV; + } else { + assert(RegionAttributeHelpers::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 (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell; + } else if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell; + } else { + assert(RegionAttributeHelpers::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::RegionId RegionId; + + /** + * Average pressure + * + */ + double + pressure(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 rmap_; + + /** + * Derived property attributes for each active region. + */ + struct Attributes { + Attributes() + : pressure(0.0) + , pv(0.0) + + {} + + double pressure; + double pv; + + }; + + RegionAttributeHelpers::RegionAttributes attr_; + + }; + } // namespace RegionAverageCalculator +} // namespace Opm + +#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */ diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 59b8c817d..6393cf435 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -211,7 +211,7 @@ namespace WellGroupHelpers return; const auto [name, number] = *gpm->region(); - const double error = gpm->pressure_target() - regional_values.fpr(number); + const double error = gpm->pressure_target() - regional_values.pressure(number); double current_rate = 0.0; const auto& pu = well_state.phaseUsage(); double sign = 1.0;