From 9d2f26f7e88081027f3f160cfc83c0404dac51fc Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 1 Sep 2021 15:04:54 +0200 Subject: [PATCH 1/4] Add support for gpmaint --- opm/simulators/wells/BlackoilWellModel.hpp | 6 + .../wells/BlackoilWellModelGeneric.cpp | 28 ++- .../wells/BlackoilWellModel_impl.hpp | 14 +- opm/simulators/wells/GroupState.cpp | 39 +++ opm/simulators/wells/GroupState.hpp | 11 + opm/simulators/wells/RateConverter.hpp | 233 ++++++++++++++++++ opm/simulators/wells/TargetCalculator.cpp | 24 +- opm/simulators/wells/TargetCalculator.hpp | 10 +- opm/simulators/wells/WellGroupHelpers.cpp | 43 +++- opm/simulators/wells/WellGroupHelpers.hpp | 90 +++++++ opm/simulators/wells/WellInterfaceEval.cpp | 4 +- .../wells/WellInterfaceFluidSystem.cpp | 4 +- 12 files changed, 487 insertions(+), 19 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index f1e315827..e685ceebe 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -132,6 +132,10 @@ namespace Opm { using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; + // For computing average pressured used by gpmaint + using AverageRegionalPressureType = RateConverter:: + AverageRegionalPressure >; + BlackoilWellModel(Simulator& ebosSimulator); void init(); @@ -317,6 +321,8 @@ namespace Opm { bool alternative_well_rate_init_{}; std::unique_ptr rateConverter_{}; + std::unique_ptr gpmaint_{}; + SimulatorReportSingle last_report_{}; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 59af13092..848c2cdfc 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -450,10 +450,15 @@ checkGroupInjectionConstraints(const Group& group, current_rate = comm_.sum(current_rate); const auto& controls = group.injectionControls(phase, this->summaryState_); - if (controls.surface_max_rate < current_rate) { + double target = controls.surface_max_rate; + + if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE) && this->groupState().has_gpmaint_target(group.name())) + target = this->groupState().gpmaint_target(group.name()); + + if (target < current_rate) { double scale = 1.0; if (current_rate > 1e-12) - scale = controls.surface_max_rate / current_rate; + scale = target / current_rate; return std::make_pair(Group::InjectionCMode::RATE, scale); } } @@ -468,10 +473,15 @@ checkGroupInjectionConstraints(const Group& group, current_rate = comm_.sum(current_rate); const auto& controls = group.injectionControls(phase, this->summaryState_); - if (controls.resv_max_rate < current_rate) { + double target = controls.resv_max_rate; + + if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV) && this->groupState().has_gpmaint_target(group.name())) + target = this->groupState().gpmaint_target(group.name()); + + if (target < current_rate) { double scale = 1.0; if (current_rate > 1e-12) - scale = controls.resv_max_rate / current_rate; + scale = target / current_rate; return std::make_pair(Group::InjectionCMode::RESV, scale); } } @@ -637,10 +647,14 @@ checkGroupProductionConstraints(const Group& group, // sum over all nodes current_rate = comm_.sum(current_rate); - if (controls.resv_target < current_rate ) { + double target = controls.resv_target; + if (group.has_gpmaint_control(Group::ProductionCMode::RESV) && this->groupState().has_gpmaint_target(group.name())) + target = this->groupState().gpmaint_target(group.name()); + + if ( target < current_rate ) { double scale = 1.0; if (current_rate > 1e-12) - scale = controls.resv_target / current_rate; + scale = target / current_rate; return std::make_pair(Group::ProductionCMode::RESV, scale); } } @@ -1521,6 +1535,8 @@ updateAndCommunicateGroupData(const int reportStepIdx, WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + WellGroupHelpers::updateSurfaceRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); // We use the rates from the previous time-step to reduce oscillations diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 6652363ff..9d9d5326e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -207,6 +207,13 @@ namespace Opm { std::vector(local_num_cells_, 0))); rateConverter_->template defineState(ebosSimulator_); + // Compute reservoir volumes for RESV controls. + const auto& fp = this->eclState_.fieldProps(); + const auto fipnum = fp.get_int("FIPNUM"); + gpmaint_.reset(new AverageRegionalPressureType (phase_usage_, + fipnum)); + gpmaint_->template defineState(ebosSimulator_); + { const auto& sched_state = this->schedule()[timeStepIdx]; // update VFP properties @@ -458,6 +465,12 @@ namespace Opm { // update the rate converter with current averages pressures etc in rateConverter_->template defineState(ebosSimulator_); + gpmaint_->template defineState(ebosSimulator_); + const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx); + WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup, + schedule_, *gpmaint_, reportStepIdx, dt, this->wellState(), this->groupState()); + + // calculate the well potentials try { updateWellPotentials(reportStepIdx, @@ -472,7 +485,6 @@ namespace Opm { updateWellTestState(simulationTime, wellTestState_); // check group sales limits at the end of the timestep - const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); checkGconsaleLimits(fieldGroup, this->wellState(), ebosSimulator_.episodeIndex(), local_deferredLogger); diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index f663f5e99..ae7be4672 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -107,6 +107,27 @@ const std::vector& GroupState::injection_reduction_rates(const std::stri return group_iter->second; } +//------------------------------------------------------------------------- + +bool GroupState::has_injection_surface_rates(const std::string& gname) const { + auto group_iter = this->inj_surface_rates.find(gname); + return (group_iter != this->inj_surface_rates.end()); +} + +void GroupState::update_injection_surface_rates(const std::string& gname, const std::vector& rates) { + if (rates.size() != this->num_phases) + throw std::logic_error("Wrong number of phases"); + + this->inj_surface_rates[gname] = rates; +} + +const std::vector& GroupState::injection_surface_rates(const std::string& gname) const { + auto group_iter = this->inj_surface_rates.find(gname); + if (group_iter == this->inj_surface_rates.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} //------------------------------------------------------------------------- @@ -246,6 +267,24 @@ GPMaint::State& GroupState::gpmaint(const std::string& gname) { } +//------------------------------------------------------------------------- + +void GroupState::update_gpmaint_target(const std::string& gname, double target) { + this->m_gpmaint_target[gname] = target; +} + +double GroupState::gpmaint_target(const std::string& gname) const { + auto group_iter = this->m_gpmaint_target.find(gname); + if (group_iter == this->m_gpmaint_target.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} + +bool GroupState::has_gpmaint_target(const std::string& gname) const { + return (this->m_gpmaint_target.count(gname) > 0); +} + //------------------------------------------------------------------------- namespace { diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index e8504963b..2c53e618b 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -52,6 +52,11 @@ public: void update_injection_reservoir_rates(const std::string& gname, const std::vector& rates); const std::vector& injection_reservoir_rates(const std::string& gname) const; + bool has_injection_surface_rates(const std::string& gname) const; + void update_injection_surface_rates(const std::string& gname, const std::vector& rates); + const std::vector& injection_surface_rates(const std::string& gname) const; + + void update_injection_rein_rates(const std::string& gname, const std::vector& rates); const std::vector& injection_rein_rates(const std::string& gname) const; @@ -65,6 +70,10 @@ public: double grat_sales_target(const std::string& gname) const; bool has_grat_sales_target(const std::string& gname) const; + void update_gpmaint_target(const std::string& gname, double target); + double gpmaint_target(const std::string& gname) const; + bool has_gpmaint_target(const std::string& gname) const; + bool has_production_control(const std::string& gname) const; void production_control(const std::string& gname, Group::ProductionCMode cmode); Group::ProductionCMode production_control(const std::string& gname) const; @@ -161,11 +170,13 @@ private: std::map production_controls; std::map> prod_red_rates; std::map> inj_red_rates; + std::map> inj_surface_rates; std::map> inj_resv_rates; std::map> inj_potentials; std::map> inj_rein_rates; std::map inj_vrep_rate; std::map m_grat_sales_target; + std::map m_gpmaint_target; std::map, Group::InjectionCMode> injection_controls; diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 01de81bcd..ab25d164d 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -197,6 +197,32 @@ namespace Opm { 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. * @@ -855,6 +881,213 @@ namespace Opm { 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_; + + }; } // namespace RateConverter } // namespace Opm diff --git a/opm/simulators/wells/TargetCalculator.cpp b/opm/simulators/wells/TargetCalculator.cpp index 7a5904cf0..ac52d88c5 100644 --- a/opm/simulators/wells/TargetCalculator.cpp +++ b/opm/simulators/wells/TargetCalculator.cpp @@ -40,13 +40,19 @@ namespace WellGroupHelpers TargetCalculator::TargetCalculator(const Group::ProductionCMode cmode, const PhaseUsage& pu, const std::vector& resv_coeff, - const double group_grat_target_from_sales) + const double group_grat_target_from_sales, + const std::string& group_name, + const GroupState& group_state, + const bool use_gpmaint) : cmode_(cmode) , pu_(pu) , resv_coeff_(resv_coeff) , group_grat_target_from_sales_(group_grat_target_from_sales) -{ + , group_name_(group_name) + , group_state_(group_state) + , use_gpmaint_(use_gpmaint) +{ } template @@ -107,7 +113,12 @@ double TargetCalculator::groupTarget(const Group::ProductionControls ctrl) const case Group::ProductionCMode::LRAT: return ctrl.liquid_target; case Group::ProductionCMode::RESV: + { + if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_)) + return this->group_state_.gpmaint_target(this->group_name_); + return ctrl.resv_target; + } default: // Should never be here. assert(false); @@ -142,6 +153,7 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode const double sales_target, const GroupState& group_state, const Phase& injection_phase, + const bool use_gpmaint, DeferredLogger& deferred_logger) : cmode_(cmode) , pu_(pu) @@ -149,6 +161,8 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode , group_name_(group_name) , sales_target_(sales_target) , group_state_(group_state) + , use_gpmaint_(use_gpmaint) + { // initialize to avoid warning pos_ = pu.phase_pos[BlackoilPhases::Aqua]; @@ -182,8 +196,14 @@ double InjectionTargetCalculator::groupTarget(const Group::InjectionControls& ct { switch (cmode_) { case Group::InjectionCMode::RATE: + if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_)) + return this->group_state_.gpmaint_target(this->group_name_); + return ctrl.surface_max_rate; case Group::InjectionCMode::RESV: + if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_)) + return this->group_state_.gpmaint_target(this->group_name_) / resv_coeff_[pos_]; + return ctrl.resv_max_rate / resv_coeff_[pos_]; case Group::InjectionCMode::REIN: { double production_rate = this->group_state_.injection_rein_rates(ctrl.reinj_group)[pos_]; diff --git a/opm/simulators/wells/TargetCalculator.hpp b/opm/simulators/wells/TargetCalculator.hpp index fe2566759..40b8bbd3c 100644 --- a/opm/simulators/wells/TargetCalculator.hpp +++ b/opm/simulators/wells/TargetCalculator.hpp @@ -45,7 +45,10 @@ namespace WellGroupHelpers TargetCalculator(const Group::ProductionCMode cmode, const PhaseUsage& pu, const std::vector& resv_coeff, - const double group_grat_target_from_sales); + const double group_grat_target_from_sales, + const std::string& group_name, + const GroupState& group_state, + const bool use_gpmaint); template RateType calcModeRateFromRates(const std::vector& rates) const @@ -65,6 +68,9 @@ namespace WellGroupHelpers const PhaseUsage& pu_; const std::vector& resv_coeff_; const double group_grat_target_from_sales_; + const std::string& group_name_; + const GroupState& group_state_; + bool use_gpmaint_; }; /// Based on a group control mode, extract or calculate rates, and @@ -79,6 +85,7 @@ namespace WellGroupHelpers const double sales_target, const GroupState& group_state, const Phase& injection_phase, + const bool use_gpmaint, DeferredLogger& deferred_logger); template @@ -98,6 +105,7 @@ namespace WellGroupHelpers const std::string& group_name_; double sales_target_; const GroupState& group_state_; + bool use_gpmaint_; int pos_; GuideRateModel::Target target_; }; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index b40c7957c..b4133b6e0 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -21,13 +21,10 @@ #include #include -#include #include #include -#include #include #include -#include #include #include #include @@ -163,6 +160,17 @@ namespace WellGroupHelpers group_state.production_control(group.name(), controls.cmode); } + if (group.has_gpmaint_control(Group::ProductionCMode::RESV)) { + group_state.production_control(group.name(), Group::ProductionCMode::RESV); + } + for (Phase phase : all) { + if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE)) { + group_state.injection_control(group.name(), phase, Group::InjectionCMode::RATE); + } else if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV)) { + group_state.injection_control(group.name(), phase, Group::InjectionCMode::RESV); + } + } + if (schedule[reportStepIdx].gconsale().has(group.name())) { group_state.injection_control(group.name(), Phase::GAS, Group::InjectionCMode::SALE); } @@ -546,6 +554,31 @@ namespace WellGroupHelpers group_state.update_injection_reservoir_rates(group.name(), resv); } + void updateSurfaceRatesInjectionGroups(const Group& group, + const Schedule& schedule, + const int reportStepIdx, + const WellState& wellStateNupcol, + WellState& wellState, + GroupState& group_state) + { + for (const std::string& groupName : group.groups()) { + const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); + updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellStateNupcol, wellState, group_state); + } + const int np = wellState.numPhases(); + std::vector rates(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + rates[phase] = sumWellPhaseRates(false, + group, + schedule, + wellStateNupcol, + reportStepIdx, + phase, + /*isInjector*/ true); + } + group_state.update_injection_surface_rates(group.name(), rates); + } + void updateWellRates(const Group& group, const Schedule& schedule, const int reportStepIdx, @@ -1118,7 +1151,7 @@ namespace WellGroupHelpers if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); - TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales); + TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl)); FractionCalculator fcalc(schedule, wellState, group_state, reportStepIdx, guideRate, tcalc.guideTargetMode(), pu, true, Phase::OIL); auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); }; @@ -1245,7 +1278,7 @@ namespace WellGroupHelpers const auto& gconsale = schedule[reportStepIdx].gconsale().get(group.name(), summaryState); sales_target = gconsale.sales_target; } - InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger); + InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger); FractionCalculator fcalc(schedule, wellState, group_state, reportStepIdx, guideRate, tcalc.guideTargetMode(), pu, false, injectionPhase); auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); }; diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 2aea55564..59b8c817d 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -22,6 +22,11 @@ #define OPM_WELLGROUPHELPERS_HEADER_INCLUDED #include +#include +#include +#include +#include +#include #include #include @@ -151,6 +156,13 @@ namespace WellGroupHelpers WellState& wellState, GroupState& group_state); + void updateSurfaceRatesInjectionGroups(const Group& group, + const Schedule& schedule, + const int reportStepIdx, + const WellState& wellStateNupcol, + WellState& wellState, + GroupState& group_state); + void updateWellRates(const Group& group, const Schedule& schedule, const int reportStepIdx, @@ -181,6 +193,84 @@ namespace WellGroupHelpers WellState& wellState, GroupState& group_state); + template + void updateGpMaintTargetForGroups(const Group& group, + const Schedule& schedule, + const RegionalValues& regional_values, + const int reportStepIdx, + const double dt, + const WellState& well_state, + GroupState& group_state) + { + for (const std::string& groupName : group.groups()) { + const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); + updateGpMaintTargetForGroups(groupTmp, schedule, regional_values, reportStepIdx, dt, well_state, group_state); + } + const auto& gpm = group.gpmaint(); + if(!gpm) + return; + + const auto [name, number] = *gpm->region(); + const double error = gpm->pressure_target() - regional_values.fpr(number); + double current_rate = 0.0; + const auto& pu = well_state.phaseUsage(); + double sign = 1.0; + switch (gpm->flow_target()) { + case Opm::GPMaint::FlowTarget::RESV_PROD: + { + current_rate = -group_state.injection_vrep_rate(group.name()); + sign = -1.0; + break; + } + case Opm::GPMaint::FlowTarget::RESV_OINJ: + { + if (pu.phase_used[BlackoilPhases::Liquid]) + current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Liquid]]; + + break; + } + case Opm::GPMaint::FlowTarget::RESV_WINJ: + { + if (pu.phase_used[BlackoilPhases::Aqua]) + current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Aqua]]; + + break; + } + case Opm::GPMaint::FlowTarget::RESV_GINJ: + { + if (pu.phase_used[BlackoilPhases::Vapour]) + current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Vapour]]; + break; + } + case Opm::GPMaint::FlowTarget::SURF_OINJ: + { + if (pu.phase_used[BlackoilPhases::Liquid]) + current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Liquid]]; + + break; + } + case Opm::GPMaint::FlowTarget::SURF_WINJ: + { + if (pu.phase_used[BlackoilPhases::Aqua]) + current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Aqua]]; + + break; + } + case Opm::GPMaint::FlowTarget::SURF_GINJ: + { + if (pu.phase_used[BlackoilPhases::Vapour]) + current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Vapour]]; + + break; + } + default: + throw std::invalid_argument("Invalid Flow target type in GPMAINT"); + } + + auto& gpmaint_state = group_state.gpmaint(group.name()); + double rate = gpm->rate(gpmaint_state, current_rate, error, dt); + group_state.update_gpmaint_target(group.name(), std::max(0.0, sign * rate)); + } std::map computeNetworkPressures(const Opm::Network::ExtNetwork& network, const WellState& well_state, diff --git a/opm/simulators/wells/WellInterfaceEval.cpp b/opm/simulators/wells/WellInterfaceEval.cpp index 72b89e2f7..037d83862 100644 --- a/opm/simulators/wells/WellInterfaceEval.cpp +++ b/opm/simulators/wells/WellInterfaceEval.cpp @@ -136,7 +136,7 @@ getGroupInjectionControl(const Group& group, const auto& gconsale = schedule[baseif_.currentStep()].gconsale().get(group.name(), summaryState); sales_target = gconsale.sales_target; } - WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger); + WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase); auto localFraction = [&](const std::string& child) { @@ -230,7 +230,7 @@ getGroupProductionControl(const Group& group, if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); - WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales); + WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl)); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL); auto localFraction = [&](const std::string& child) { diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index 7e70b3bee..ce7e24d2d 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -1002,7 +1002,7 @@ getGroupInjectionTargetRate(const Group& group, const auto& gconsale = schedule[currentStep()].gconsale().get(group.name(), summaryState); sales_target = gconsale.sales_target; } - WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger); + WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, currentStep(), guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase); auto localFraction = [&](const std::string& child) { @@ -1072,7 +1072,7 @@ getGroupProductionTargetRate(const Group& group, if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); - WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales); + WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl)); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, currentStep(), guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL); auto localFraction = [&](const std::string& child) { From f9832d8830a7a90a8a85cfd0888d91bfccc3a42f Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 6 Sep 2021 14:20:40 +0200 Subject: [PATCH 2/4] split RateConverter --- CMakeLists_files.cmake | 2 + opm/simulators/wells/BlackoilWellModel.hpp | 3 +- opm/simulators/wells/RateConverter.hpp | 636 +----------------- .../wells/RegionAttributeHelpers.hpp | 409 +++++++++++ .../wells/RegionAverageCalculator.hpp | 258 +++++++ opm/simulators/wells/WellGroupHelpers.hpp | 2 +- 6 files changed, 704 insertions(+), 606 deletions(-) create mode 100644 opm/simulators/wells/RegionAttributeHelpers.hpp create mode 100644 opm/simulators/wells/RegionAverageCalculator.hpp 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; From f497fa0fd7e4578e6846ff70fe90950d49086668 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 6 Sep 2021 14:54:26 +0200 Subject: [PATCH 3/4] add support for region = 0 i.e. field --- .../wells/RegionAttributeHelpers.hpp | 51 ++++++++++++------- .../wells/RegionAverageCalculator.hpp | 13 +++++ 2 files changed, 45 insertions(+), 19 deletions(-) diff --git a/opm/simulators/wells/RegionAttributeHelpers.hpp b/opm/simulators/wells/RegionAttributeHelpers.hpp index e65e7a32e..e8a825ac4 100644 --- a/opm/simulators/wells/RegionAttributeHelpers.hpp +++ b/opm/simulators/wells/RegionAttributeHelpers.hpp @@ -140,6 +140,27 @@ namespace Opm { typename Select::RegionIDParameter ::value>::type; + using ID = + typename std::remove_reference::type; + + /** + * 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 AttributeMap = + std::unordered_map>; + + /** * Constructor. * @@ -213,6 +234,17 @@ namespace Opm { } } + /** + * Request read-only access to region's attributes. + * + * + * \return Read-only access to all regions attributes. + */ + const AttributeMap& attributes() const + { + return attr_; + } + /** * Request read-only access to region's attributes. @@ -241,25 +273,6 @@ namespace Opm { } 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_; diff --git a/opm/simulators/wells/RegionAverageCalculator.hpp b/opm/simulators/wells/RegionAverageCalculator.hpp index df7663944..14356232f 100644 --- a/opm/simulators/wells/RegionAverageCalculator.hpp +++ b/opm/simulators/wells/RegionAverageCalculator.hpp @@ -218,6 +218,19 @@ namespace Opm { double pressure(const RegionId r) const { + if (r == 0 ) // region 0 is the whole field + { + double pressure = 0.0; + int num_active_regions = 0; + for (const auto& attr : attr_.attributes()) { + const auto& value = *attr.second; + const auto& ra = value.attr_; + pressure += ra.pressure; + num_active_regions ++; + } + return pressure / num_active_regions; + } + const auto& ra = attr_.attributes(r); return ra.pressure; } From adda160098065062255369d7d100ab9de83537a7 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 6 Sep 2021 15:02:20 +0200 Subject: [PATCH 4/4] clean-up --- opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModelGeneric.cpp | 2 +- .../wells/BlackoilWellModel_impl.hpp | 16 +++++++------- .../wells/RegionAttributeHelpers.hpp | 20 ++++-------------- .../wells/RegionAverageCalculator.hpp | 21 +++++-------------- opm/simulators/wells/WellGroupHelpers.cpp | 13 ++++++------ opm/simulators/wells/WellGroupHelpers.hpp | 9 ++++---- 7 files changed, 28 insertions(+), 55 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 83bab1b1a..7399c1de8 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -322,7 +322,7 @@ namespace Opm { bool alternative_well_rate_init_{}; std::unique_ptr rateConverter_{}; - std::unique_ptr gpmaint_{}; + std::unique_ptr regionalAveragePressureCalculator_{}; SimulatorReportSingle last_report_{}; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 848c2cdfc..c5c34a49e 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1535,7 +1535,7 @@ updateAndCommunicateGroupData(const int reportStepIdx, WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); - WellGroupHelpers::updateSurfaceRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + WellGroupHelpers::updateSurfaceRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, this->groupState()); WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 9d9d5326e..e420595ac 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -203,16 +203,14 @@ namespace Opm { WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState()); // Compute reservoir volumes for RESV controls. - rateConverter_.reset(new RateConverterType (phase_usage_, - std::vector(local_num_cells_, 0))); + rateConverter_.reset(new RateConverterType (phase_usage_, std::vector(local_num_cells_, 0))); rateConverter_->template defineState(ebosSimulator_); - // Compute reservoir volumes for RESV controls. + // Compute regional average pressures used by gpmaint const auto& fp = this->eclState_.fieldProps(); - const auto fipnum = fp.get_int("FIPNUM"); - gpmaint_.reset(new AverageRegionalPressureType (phase_usage_, - fipnum)); - gpmaint_->template defineState(ebosSimulator_); + const auto& fipnum = fp.get_int("FIPNUM"); + regionalAveragePressureCalculator_.reset(new AverageRegionalPressureType (phase_usage_,fipnum)); + regionalAveragePressureCalculator_->template defineState(ebosSimulator_); { const auto& sched_state = this->schedule()[timeStepIdx]; @@ -465,10 +463,10 @@ namespace Opm { // update the rate converter with current averages pressures etc in rateConverter_->template defineState(ebosSimulator_); - gpmaint_->template defineState(ebosSimulator_); + regionalAveragePressureCalculator_->template defineState(ebosSimulator_); const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx); WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup, - schedule_, *gpmaint_, reportStepIdx, dt, this->wellState(), this->groupState()); + schedule_, *regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState()); // calculate the well potentials diff --git a/opm/simulators/wells/RegionAttributeHelpers.hpp b/opm/simulators/wells/RegionAttributeHelpers.hpp index e8a825ac4..b7d2d67a5 100644 --- a/opm/simulators/wells/RegionAttributeHelpers.hpp +++ b/opm/simulators/wells/RegionAttributeHelpers.hpp @@ -211,26 +211,14 @@ namespace Opm { bool has(const RegionID reg) const { - const auto& i = attr_.find(reg); - - if (i == attr_.end()) { - return false; - } - - return true; + return this->attr_.find(reg) != this->attr_.end(); } 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; + auto [pos, inserted] = this->attr_.try_emplace(r, std::make_unique(attr)); + if (inserted) { + pos->second->cell_ = -1; // NOT -1.0 -- "cell_" is 'int' } } diff --git a/opm/simulators/wells/RegionAverageCalculator.hpp b/opm/simulators/wells/RegionAverageCalculator.hpp index 14356232f..430811bbb 100644 --- a/opm/simulators/wells/RegionAverageCalculator.hpp +++ b/opm/simulators/wells/RegionAverageCalculator.hpp @@ -68,7 +68,7 @@ namespace Opm { * deck. */ AverageRegionalPressure(const PhaseUsage& phaseUsage, - const Region& region) + const Region& region) : phaseUsage_(phaseUsage) , rmap_ (region) , attr_ (rmap_, Attributes()) @@ -91,13 +91,13 @@ namespace Opm { numRegions = std::max(numRegions, reg); } numRegions = comm.max(numRegions); - for (int reg = 0; reg < numRegions; ++ reg) { + for (int reg = 1; 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) { + for (int reg = 1; reg <= numRegions ; ++ reg) { auto& ra = attr_.attributes(reg); ra.pressure = 0.0; ra.pv = 0.0; @@ -110,22 +110,11 @@ namespace Opm { // quantities for hydrocarbon volume average std::unordered_map attributes_hpv; - for (int reg = 0; reg < numRegions; ++ reg) { + for (int reg = 1; 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(); @@ -184,7 +173,7 @@ namespace Opm { } } - for (int reg = 0; reg < numRegions; ++ reg) { + for (int reg = 1; 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? diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index b4133b6e0..c935dee14 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -555,17 +555,16 @@ namespace WellGroupHelpers } void updateSurfaceRatesInjectionGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellStateNupcol, - WellState& wellState, - GroupState& group_state) + const Schedule& schedule, + const int reportStepIdx, + const WellState& wellStateNupcol, + GroupState& group_state) { for (const std::string& groupName : group.groups()) { const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellStateNupcol, wellState, group_state); + updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellStateNupcol, group_state); } - const int np = wellState.numPhases(); + const int np = wellStateNupcol.numPhases(); std::vector rates(np, 0.0); for (int phase = 0; phase < np; ++phase) { rates[phase] = sumWellPhaseRates(false, diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 6393cf435..82e60d2f7 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -157,11 +157,10 @@ namespace WellGroupHelpers GroupState& group_state); void updateSurfaceRatesInjectionGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellStateNupcol, - WellState& wellState, - GroupState& group_state); + const Schedule& schedule, + const int reportStepIdx, + const WellState& wellStateNupcol, + GroupState& group_state); void updateWellRates(const Group& group, const Schedule& schedule,