diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 380613f59..3051c8d0e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -218,7 +218,6 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/PerforationData.hpp opm/simulators/wells/RateConverter.hpp opm/simulators/utils/readDeck.hpp - opm/simulators/wells/SimFIBODetails.hpp opm/simulators/wells/TargetCalculator.hpp opm/simulators/wells/WellConnectionAuxiliaryModule.hpp opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp diff --git a/examples/printvfp.cpp b/examples/printvfp.cpp index eb8ce71a5..e2b357235 100644 --- a/examples/printvfp.cpp +++ b/examples/printvfp.cpp @@ -93,10 +93,11 @@ int main(int argc, char** argv) } Setup setup(argv[1]); -// const int table_id = 1; +// const int table_id = 1; const int table_id = 4; const double wct = 0.0; - const double gor = 35.2743; +// const double gor = 35.2743; + const double gor = 0.0; const double alq = 0.0; const int n = 51; const double m3pd = unit::cubic(unit::meter)/unit::day; diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index f776617b4..27d304cd6 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -145,6 +145,10 @@ template struct MaxInnerIterWells { using type = UndefinedProperty; }; +template +struct AlternativeWellRateInit { + using type = UndefinedProperty; +}; template struct DbhpMaxRel { @@ -251,6 +255,10 @@ struct MaxInnerIterWells { static constexpr int value = 50; }; template +struct AlternativeWellRateInit { + static constexpr bool value = false; +}; +template struct StrictInnerIterMsWells { static constexpr int value = 40; }; @@ -432,6 +440,7 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, int, StrictInnerIterMsWells, "Number of inner iterations for multi-segment wells with strict tolerance"); EWOMS_REGISTER_PARAM(TypeTag, bool, UseInnerIterationsWells, "Use nested iterations for standard wells"); EWOMS_REGISTER_PARAM(TypeTag, int, MaxInnerIterWells, "Maximum number of inner iterations for standard wells"); + EWOMS_REGISTER_PARAM(TypeTag, bool, AlternativeWellRateInit, "Use alternative well rate initialization procedure"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, RegularizationFactorMsw, "Regularization factor for ms wells"); EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxSinglePrecisionDays, "Maximum time step size where single precision floating point arithmetic can be used solving for the linear systems of equations"); EWOMS_REGISTER_PARAM(TypeTag, int, MaxStrictIter, "Maximum number of Newton iterations before relaxed tolerances are used for the CNV convergence criterion"); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 9429febab..e638e1002 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -196,8 +196,8 @@ namespace Opm { { auto grp_nwrk_values = ::Opm::data::GroupAndNetworkValues{}; - this->assignGroupValues(reportStepIdx, sched, - grp_nwrk_values.groupData); + this->assignGroupValues(reportStepIdx, sched, grp_nwrk_values.groupData); + this->assignNodeValues(reportStepIdx, sched, grp_nwrk_values.nodeData); return grp_nwrk_values; } @@ -307,6 +307,7 @@ namespace Opm { bool initial_step_; bool report_step_starts_; bool glift_debug = false; + bool alternative_well_rate_init_; std::unique_ptr rateConverter_; std::unique_ptr> vfp_properties_; @@ -315,6 +316,8 @@ namespace Opm { WellTestState wellTestState_; std::unique_ptr guideRate_; + std::map node_pressures_; // Storing network pressures for output. + // used to better efficiency of calcuation mutable BVector scaleAddRes_; @@ -355,6 +358,7 @@ namespace Opm { void updateWellControls(Opm::DeferredLogger& deferred_logger, const bool checkGroupControls); void updateAndCommunicateGroupData(); + void updateNetworkPressures(); // setting the well_solutions_ based on well_state. void updatePrimaryVariables(Opm::DeferredLogger& deferred_logger); @@ -452,6 +456,10 @@ namespace Opm { const Schedule& sched, std::map& gvalues) const; + void assignNodeValues(const int reportStepIdx, + const Schedule& sched, + std::map& gvalues) const; + std::unordered_map calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 6d368c680..f808f7ef6 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -21,7 +21,6 @@ */ #include -#include #include #include @@ -64,6 +63,8 @@ namespace Opm { value); return candidate == parallel_wells.end() || *candidate != value; }; + + alternative_well_rate_init_ = EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit); } template @@ -329,6 +330,8 @@ namespace Opm { BlackoilWellModel:: beginTimeStep() { + updatePerforationIntensiveQuantities(); + Opm::DeferredLogger local_deferredLogger; well_state_ = previous_well_state_; @@ -392,6 +395,13 @@ namespace Opm { local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg); } + if (alternative_well_rate_init_) { + // Update the well rates to match state, if only single-phase rates. + for (auto& well : well_container_) { + well->updateWellStateRates(ebosSimulator_, well_state_, local_deferredLogger); + } + } + //compute well guideRates const auto& comm = ebosSimulator_.vanguard().grid().comm(); WellGroupHelpers::updateGuideRatesForWells(schedule(), phase_usage_, reportStepIdx, simulationTime, well_state_, comm, guideRate_.get()); @@ -1219,6 +1229,8 @@ namespace Opm { updateAndCommunicateGroupData(); + updateNetworkPressures(); + std::set switched_wells; std::set switched_groups; @@ -1257,6 +1269,37 @@ namespace Opm { + template + void + BlackoilWellModel:: + updateNetworkPressures() + { + // Get the network and return if inactive. + const int reportStepIdx = ebosSimulator_.episodeIndex(); + const auto& network = schedule().network(reportStepIdx); + if (!network.active()) { + return; + } + node_pressures_ = WellGroupHelpers::computeNetworkPressures(network, well_state_, *(vfp_properties_->getProd())); + + // Set the thp limits of wells + for (auto& well : well_container_) { + // Producers only, since we so far only support the + // "extended" network model (properties defined by + // BRANPROP and NODEPROP) which only applies to producers. + if (well->isProducer()) { + const auto it = node_pressures_.find(well->wellEcl().groupName()); + if (it != node_pressures_.end()) { + // The well belongs to a group with has a network pressure constraint, + // set the dynamic THP constraint of the well accordingly. + well->setDynamicThpLimit(it->second); + } + } + } + } + + + template void @@ -2503,6 +2546,19 @@ namespace Opm { } } + template + void + BlackoilWellModel:: + assignNodeValues(const int reportStepIdx, + const Schedule& sched, + std::map& nodevalues) const + { + nodevalues.clear(); + for (const auto& [node, pressure] : node_pressures_) { + nodevalues.emplace(node, data::NodeData{pressure}); + } + } + template std::unordered_map BlackoilWellModel:: diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index c0694e762..e667d392b 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -176,6 +176,9 @@ namespace Opm int numberOfPerforations() const; + virtual std::vector computeCurrentWellRates(const Simulator& ebosSimulator, + DeferredLogger& deferred_logger) const override; + protected: int number_segments_; @@ -333,6 +336,7 @@ namespace Opm void computePerfRatePressure(const IntensiveQuantities& int_quants, const std::vector& mob_perfcells, + const double Tw, const int seg, const int perf, const EvalWell& segment_pressure, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 9fd5b5bcc..c78ff0174 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1321,6 +1321,7 @@ namespace Opm MultisegmentWell:: computePerfRatePressure(const IntensiveQuantities& int_quants, const std::vector& mob_perfcells, + const double Tw, const int seg, const int perf, const EvalWell& segment_pressure, @@ -1376,7 +1377,7 @@ namespace Opm // compute component volumetric rates at standard conditions for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell cq_p = - well_index_[perf] * (mob_perfcells[comp_idx] * drawdown); + const EvalWell cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown); cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p; } @@ -1401,7 +1402,7 @@ namespace Opm } // injection perforations total volume rates - const EvalWell cqt_i = - well_index_[perf] * (total_mob * drawdown); + const EvalWell cqt_i = - Tw * (total_mob * drawdown); // compute volume ratio between connection and at standard conditions EvalWell volume_ratio = 0.0; @@ -2026,13 +2027,13 @@ namespace Opm const auto& controls = well.injectionControls(summaryState); const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit) - dp; + return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState)) - dp; } else if (well.isProducer()) { const auto& controls = well.productionControls(summaryState); const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit, controls.alq_value) - dp; + return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), controls.alq_value) - dp; } else { OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger); @@ -2602,12 +2603,13 @@ namespace Opm const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob); + const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); + const double Tw = well_index_[perf] * trans_mult; std::vector cq_s(num_components_, 0.0); EvalWell perf_press; double perf_dis_gas_rate = 0.; double perf_vap_oil_rate = 0.; - - computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + computePerfRatePressure(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); // updating the solution gas rate and solution oil rate if (this->isProducer()) { @@ -3306,11 +3308,12 @@ namespace Opm const auto& table = *(vfp_properties_->getProd()->getTable(controls.vfp_table_number)); const double vfp_ref_depth = table.getDatumDepth(); const double rho = segment_densities_[0].value(); // Use the density at the top perforation. + const double thp_limit = this->getTHPConstraint(summary_state); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - auto fbhp = [this, &controls, dp](const std::vector& rates) { + auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { assert(rates.size() == 3); return this->vfp_properties_->getProd() - ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit, controls.alq_value) - dp; + ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, controls.alq_value) - dp; }; // Make the flo() function. @@ -3529,11 +3532,12 @@ namespace Opm const auto& table = *(vfp_properties_->getInj()->getTable(controls.vfp_table_number)); const double vfp_ref_depth = table.getDatumDepth(); const double rho = segment_densities_[0].value(); // Use the density at the top perforation. + const double thp_limit = this->getTHPConstraint(summary_state); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - auto fbhp = [this, &controls, dp](const std::vector& rates) { + auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { assert(rates.size() == 3); return this->vfp_properties_->getInj() - ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit) - dp; + ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp; }; // Make the flo() function. @@ -3817,4 +3821,47 @@ namespace Opm const double sign = mass_rate <= 0. ? 1.0 : -1.0; return sign * (friction_pressure_loss + constriction_pressure_loss); } + + + + + template + std::vector + MultisegmentWell:: + computeCurrentWellRates(const Simulator& ebosSimulator, + DeferredLogger& deferred_logger) const + { + // Calculate the rates that follow from the current primary variables. + std::vector well_q_s(num_components_, 0.0); + const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); + const int nseg = numberOfSegments(); + for (int seg = 0; seg < nseg; ++seg) { + // calculating the perforation rate for each perforation that belongs to this segment + const EvalWell seg_pressure = getSegmentPressure(seg); + for (const int perf : segment_perforations_[seg]) { + const int cell_idx = well_cells_[perf]; + const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + std::vector mob(num_components_, 0.0); + getMobility(ebosSimulator, perf, mob); + const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); + const double Tw = well_index_[perf] * trans_mult; + std::vector cq_s(num_components_, 0.0); + EvalWell perf_press; + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + computePerfRatePressure(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + for (int comp = 0; comp < num_components_; ++comp) { + well_q_s[comp] += cq_s[comp]; + } + } + } + std::vector well_q_s_noderiv(well_q_s.size()); + for (int comp = 0; comp < num_components_; ++comp) { + well_q_s_noderiv[comp] = well_q_s[comp].value(); + } + return well_q_s_noderiv; } + + + +} // namespace Opm diff --git a/opm/simulators/wells/SimFIBODetails.hpp b/opm/simulators/wells/SimFIBODetails.hpp deleted file mode 100644 index 396c97157..000000000 --- a/opm/simulators/wells/SimFIBODetails.hpp +++ /dev/null @@ -1,83 +0,0 @@ -/* - Copyright 2013 SINTEF ICT, Applied Mathematics. - Copyright 2014-2016 IRIS AS - Copyright 2015 Andreas Lauser - - 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_SIM_FIBO_DETAILS_HPP -#define OPM_SIM_FIBO_DETAILS_HPP - -#include -#include -#include - -#include -#include - -namespace Opm -{ - namespace SimFIBODetails { - typedef std::unordered_map WellMap; - - inline WellMap - mapWells(const std::vector< Well >& wells) - { - WellMap wmap; - - for (const auto& w : wells) - { - wmap.insert(std::make_pair(w.name(), w)); - } - - return wmap; - } - - - inline void - historyRates(const PhaseUsage& pu, - const Well::ProductionControls& p, - std::vector& rates) - { - assert (! p.prediction_mode); - assert (rates.size() == - std::vector::size_type(pu.num_phases)); - - if (pu.phase_used[ BlackoilPhases::Aqua ]) { - const std::vector::size_type - i = pu.phase_pos[ BlackoilPhases::Aqua ]; - - rates[i] = p.water_rate; - } - - if (pu.phase_used[ BlackoilPhases::Liquid ]) { - const std::vector::size_type - i = pu.phase_pos[ BlackoilPhases::Liquid ]; - - rates[i] = p.oil_rate; - } - - if (pu.phase_used[ BlackoilPhases::Vapour ]) { - const std::vector::size_type - i = pu.phase_pos[ BlackoilPhases::Vapour ]; - - rates[i] = p.gas_rate; - } - } - } // namespace SimFIBODetails -} // namespace Opm - -#endif // OPM_SIM_FIBO_DETAILS_HPP diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 85338ddf4..73383ebd2 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -296,6 +296,9 @@ namespace Opm using Base::phaseUsage; using Base::vfp_properties_; + virtual std::vector computeCurrentWellRates(const Simulator& ebosSimulator, + DeferredLogger& deferred_logger) const override; + protected: // protected functions from the Base class diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index c333c740b..e1ef46a77 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1441,7 +1441,7 @@ namespace Opm } case Well::ProducerCMode::THP: { - well_state.thp()[well_index] = controls.thp_limit; + well_state.thp()[well_index] = this->getTHPConstraint(summaryState); gliftDebug( "computing BHP from THP to update well state", deferred_logger); @@ -2900,13 +2900,13 @@ namespace Opm const auto& controls = well.injectionControls(summaryState); const double vfp_ref_depth = vfp_properties_->getInj()->getTable(controls.vfp_table_number)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit) - dp; + return vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState)) - dp; } else if (this->isProducer()) { const auto& controls = well.productionControls(summaryState); const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit, getALQ(well_state)) - dp; + return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), getALQ(well_state)) - dp; } else { OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger); @@ -3585,9 +3585,9 @@ namespace Opm std::optional StandardWell:: computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator, - const SummaryState& summary_state, - DeferredLogger& deferred_logger, - double alq_value) const + const SummaryState& summary_state, + DeferredLogger& deferred_logger, + double alq_value) const { // Given a VFP function returning bhp as a function of phase // rates and thp: @@ -3633,11 +3633,12 @@ namespace Opm const auto& table = *(vfp_properties_->getProd()->getTable(controls.vfp_table_number)); const double vfp_ref_depth = table.getDatumDepth(); const double rho = perf_densities_[0]; // Use the density at the top perforation. + const double thp_limit = this->getTHPConstraint(summary_state); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - auto fbhp = [this, &controls, dp, alq_value](const std::vector& rates) { + auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector& rates) { assert(rates.size() == 3); return this->vfp_properties_->getProd() - ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit, alq_value) - dp; + ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value) - dp; }; // Make the flo() function. @@ -3838,11 +3839,12 @@ namespace Opm const auto& table = *(vfp_properties_->getInj()->getTable(controls.vfp_table_number)); const double vfp_ref_depth = table.getDatumDepth(); const double rho = perf_densities_[0]; // Use the density at the top perforation. + const double thp_limit = this->getTHPConstraint(summary_state); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - auto fbhp = [this, &controls, dp](const std::vector& rates) { + auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { assert(rates.size() == 3); return this->vfp_properties_->getInj() - ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit) - dp; + ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp; }; // Make the flo() function. @@ -4028,4 +4030,38 @@ namespace Opm } + template + std::vector + StandardWell:: + computeCurrentWellRates(const Simulator& ebosSimulator, + DeferredLogger& deferred_logger) const + { + // Calculate the rates that follow from the current primary variables. + std::vector well_q_s(num_components_, {numWellEq_ + numEq, 0.}); + const EvalWell& bhp = getBhp(); + const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); + for (int perf = 0; perf < number_of_perforations_; ++perf) { + const int cell_idx = well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + std::vector mob(num_components_, {numWellEq_ + numEq, 0.}); + getMobility(ebosSimulator, perf, mob, deferred_logger); + std::vector cq_s(num_components_, {numWellEq_ + numEq, 0.}); + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); + const double Tw = well_index_[perf] * trans_mult; + computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf, + cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + for (int comp = 0; comp < num_components_; ++comp) { + well_q_s[comp] += cq_s[comp]; + } + } + std::vector well_q_s_noderiv(well_q_s.size()); + for (int comp = 0; comp < num_components_; ++comp) { + well_q_s_noderiv[comp] = well_q_s[comp].value(); + } + return well_q_s_noderiv; + } + + } // namespace Opm diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index dfe2759bd..19eb255fb 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -22,6 +22,7 @@ #include #include +#include #include namespace { @@ -585,6 +586,120 @@ namespace WellGroupHelpers wellState.setCurrentInjectionREINRates(group.name(), rein); } + + + + std::map + computeNetworkPressures(const Opm::Network::ExtNetwork& network, + const WellStateFullyImplicitBlackoil& well_state, + const VFPProdProperties& vfp_prod_props) + { + // TODO: Only dealing with production networks for now. + + if (!network.active()) { + return {}; + } + + // Fixed pressure nodes of the network are the roots of trees. + // Leaf nodes must correspond to groups in the group structure. + // Let us first find all leaf nodes of the network. We also + // create a vector of all nodes, ordered so that a child is + // always after its parent. + std::stack children; + std::set leaf_nodes; + std::vector root_to_child_nodes; + children.push(network.root().name()); + while (!children.empty()) { + const auto node = children.top(); + children.pop(); + root_to_child_nodes.push_back(node); + auto branches = network.downtree_branches(node); + if (branches.empty()) { + leaf_nodes.insert(node); + } + for (const auto& branch : branches) { + children.push(branch.downtree_node()); + } + } + assert(children.empty()); + + // Starting with the leaf nodes of the network, get the flow rates + // from the corresponding groups. + std::map> node_inflows; + for (const auto& node : leaf_nodes) { + node_inflows[node] = well_state.currentProductionGroupRates(node); + } + + // Accumulate in the network, towards the roots. Note that a + // root (i.e. fixed pressure node) can still be contributing + // flow towards other nodes in the network, i.e. a node is + // the root of a subtree. + auto child_to_root_nodes = root_to_child_nodes; + std::reverse(child_to_root_nodes.begin(), child_to_root_nodes.end()); + for (const auto& node : child_to_root_nodes) { + const auto upbranch = network.uptree_branch(node); + if (upbranch) { + // Add downbranch rates to upbranch. + std::vector& up = node_inflows[(*upbranch).uptree_node()]; + const std::vector& down = node_inflows[node]; + if (up.empty()) { + up = down; + } else { + assert (up.size() == down.size()); + for (size_t ii = 0; ii < up.size(); ++ii) { + up[ii] += down[ii]; + } + } + } + } + + // Going the other way (from roots to leafs), calculate the pressure + // at each node using VFP tables and rates. + std::map node_pressures; + for (const auto& node : root_to_child_nodes) { + auto press = network.node(node).terminal_pressure(); + if (press) { + node_pressures[node] = *press; + } else { + const auto upbranch = network.uptree_branch(node); + assert(upbranch); + const double up_press = node_pressures[(*upbranch).uptree_node()]; + const auto vfp_table = (*upbranch).vfp_table(); + if (vfp_table) { + // The rates are here positive, but the VFP code expects the + // convention that production rates are negative, so we must + // take a copy and flip signs. + auto rates = node_inflows[node]; + for (auto& r : rates) { r *= -1.0; } + assert(rates.size() == 3); + const double alq = 0.0; // TODO: Do not ignore ALQ + node_pressures[node] = vfp_prod_props.bhp(*vfp_table, + rates[BlackoilPhases::Aqua], + rates[BlackoilPhases::Liquid], + rates[BlackoilPhases::Vapour], + up_press, + alq); +#define EXTRA_DEBUG_NETWORK 0 +#if EXTRA_DEBUG_NETWORK + std::ostringstream oss; + oss << "parent: " << (*upbranch).uptree_node() << " child: " << node + << " rates = [ " << rates[0]*86400 << ", " << rates[1]*86400 << ", " << rates[2]*86400 << " ]" + << " p(parent) = " << up_press/1e5 << " p(child) = " << node_pressures[node]/1e5 << std::endl; + OpmLog::debug(oss.str()); +#endif + } else { + // Table number specified as 9999 in the deck, no pressure loss. + node_pressures[node] = up_press; + } + } + } + + return node_pressures; + } + + + + GuideRate::RateVector getRateVector(const WellStateFullyImplicitBlackoil& well_state, const PhaseUsage& pu, const std::string& name) { diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index d1ffece29..b319d3b67 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -26,10 +26,13 @@ #include #include #include +#include #include #include #include +#include +#include #include #include @@ -258,6 +261,11 @@ namespace WellGroupHelpers const WellStateFullyImplicitBlackoil& wellStateNupcol, WellStateFullyImplicitBlackoil& wellState); + std::map + computeNetworkPressures(const Opm::Network::ExtNetwork& network, + const WellStateFullyImplicitBlackoil& well_state, + const VFPProdProperties& vfp_prod_props); + GuideRate::RateVector getRateVector(const WellStateFullyImplicitBlackoil& well_state, const PhaseUsage& pu, const std::string& name); diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 717fd9777..5ffb9fba8 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -275,6 +275,18 @@ namespace Opm // update perforation water throughput based on solved water rate virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0; + /// Compute well rates based on current reservoir conditions and well variables. + /// Used in updateWellStateRates(). + virtual std::vector computeCurrentWellRates(const Simulator& ebosSimulator, + DeferredLogger& deferred_logger) const = 0; + + /// Modify the well_state's rates if there is only one nonzero rate. + /// If so, that rate is kept as is, but the others are set proportionally + /// to the rates returned by computeCurrentWellRates(). + void updateWellStateRates(const Simulator& ebosSimulator, + WellState& well_state, + DeferredLogger& deferred_logger) const; + void stopWell() { wellIsStopped_ = true; } @@ -288,6 +300,7 @@ namespace Opm void setWsolvent(const double wsolvent); + void setDynamicThpLimit(const double thp_limit); protected: @@ -384,6 +397,8 @@ namespace Opm double wsolvent_; + std::optional dynamic_thp_limit_; + const PhaseUsage& phaseUsage() const; int flowPhaseToEbosCompIdx( const int phaseIdx ) const; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 7ff34a93f..7c0fa6a77 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -328,6 +328,18 @@ namespace Opm + template + void + WellInterface:: + setDynamicThpLimit(const double thp_limit) + { + dynamic_thp_limit_ = thp_limit; + } + + + + + template double WellInterface:: @@ -403,6 +415,10 @@ namespace Opm WellInterface:: wellHasTHPConstraints(const SummaryState& summaryState) const { + if (dynamic_thp_limit_) { + return true; + } + if (well_ecl_.isInjector()) { const auto controls = well_ecl_.injectionControls(summaryState); if (controls.hasControl(Well::InjectorCMode::THP)) @@ -442,6 +458,9 @@ namespace Opm WellInterface:: getTHPConstraint(const SummaryState& summaryState) const { + if (dynamic_thp_limit_) { + return *dynamic_thp_limit_; + } if (well_ecl_.isInjector()) { const auto& controls = well_ecl_.injectionControls(summaryState); return controls.thp_limit; @@ -1531,7 +1550,7 @@ namespace Opm if (controls.hasControl(Well::InjectorCMode::THP) && currentControl != Well::InjectorCMode::THP) { - const auto& thp = controls.thp_limit; + const auto& thp = this->getTHPConstraint(summaryState); double current_thp = well_state.thp()[well_index]; if (thp < current_thp) { currentControl = Well::InjectorCMode::THP; @@ -1633,7 +1652,7 @@ namespace Opm if (controls.hasControl(Well::ProducerCMode::THP) && currentControl != Well::ProducerCMode::THP) { - const auto& thp = controls.thp_limit; + const auto& thp = this->getTHPConstraint(summaryState); double current_thp = well_state.thp()[well_index]; if (thp > current_thp) { currentControl = Well::ProducerCMode::THP; @@ -2268,4 +2287,47 @@ namespace Opm } + + + template + void + WellInterface:: + updateWellStateRates(const Simulator& ebosSimulator, + WellState& well_state, + DeferredLogger& deferred_logger) const + { + // Check if the rates of this well only are single-phase, do nothing + // if more than one nonzero rate. + int nonzero_rate_index = -1; + for (int p = 0; p < number_of_phases_; ++p) { + if (well_state.wellRates()[index_of_well_ * number_of_phases_ + p] != 0.0) { + if (nonzero_rate_index == -1) { + nonzero_rate_index = p; + } else { + // More than one nonzero rate. + return; + } + } + } + if (nonzero_rate_index == -1) { + // No nonzero rates. + return; + } + + // Calculate the rates that follow from the current primary variables. + std::vector well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger); + + // Set the currently-zero phase flows to be nonzero in proportion to well_q_s. + const double initial_nonzero_rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + nonzero_rate_index]; + const int comp_idx_nz = flowPhaseToEbosCompIdx(nonzero_rate_index); + for (int p = 0; p < number_of_phases_; ++p) { + if (p != nonzero_rate_index) { + const int comp_idx = flowPhaseToEbosCompIdx(p); + double& rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + p]; + rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]); + } + } + } + + } // namespace Opm