Merge pull request #2626 from atgeirr/network-thp

Network-controlled THP limit

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

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;

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -145,6 +145,10 @@ template<class TypeTag, class MyTypeTag>
 struct MaxInnerIterWells {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct AlternativeWellRateInit {
+    using type = UndefinedProperty;
+};
 
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
@@ -251,6 +255,10 @@ struct MaxInnerIterWells<TypeTag, TTag::FlowModelParameters> {
     static constexpr int value = 50;
 };
 template<class TypeTag>
+struct AlternativeWellRateInit<TypeTag, TTag::FlowModelParameters> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
 struct StrictInnerIterMsWells<TypeTag, TTag::FlowModelParameters> {
     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");

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<RateConverterType> rateConverter_;
             std::unique_ptr<VFPProperties<VFPInjProperties,VFPProdProperties>> vfp_properties_;
 
@@ -315,6 +316,8 @@ namespace Opm {
             WellTestState wellTestState_;
             std::unique_ptr<GuideRate> guideRate_;
 
+            std::map<std::string, double> 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<std::string, data::GroupData>& gvalues) const;
 
+            void assignNodeValues(const int                              reportStepIdx,
+                                  const Schedule&                        sched,
+                                  std::map<std::string, data::NodeData>& gvalues) const;
+
             std::unordered_map<std::string, data::GroupGuideRates>
             calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const;
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -21,7 +21,6 @@
 */
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
-#include <opm/simulators/wells/SimFIBODetails.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 
 #include <utility>
@@ -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<typename TypeTag>
@@ -329,6 +330,8 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     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<std::string> switched_wells;
         std::set<std::string> switched_groups;
 
@@ -1257,6 +1269,37 @@ namespace Opm {
 
 
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    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<typename TypeTag>
     void
@@ -2503,6 +2546,19 @@ namespace Opm {
         }
     }
 
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    assignNodeValues(const int                              reportStepIdx,
+                     const Schedule&                        sched,
+                     std::map<std::string, data::NodeData>& nodevalues) const
+    {
+        nodevalues.clear();
+        for (const auto& [node, pressure] : node_pressures_) {
+            nodevalues.emplace(node, data::NodeData{pressure});
+        }
+    }
+
     template<typename TypeTag>
     std::unordered_map<std::string, data::GroupGuideRates>
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/MultisegmentWell.hpp

@@ -176,6 +176,9 @@ namespace Opm
 
         int numberOfPerforations() const;
 
+        virtual std::vector<double> 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<EvalWell>& mob_perfcells,
+                                     const double Tw,
                                      const int seg,
                                      const int perf,
                                      const EvalWell& segment_pressure,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1321,6 +1321,7 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     computePerfRatePressure(const IntensiveQuantities& int_quants,
                             const std::vector<EvalWell>& 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<EvalWell> mob(num_components_, 0.0);
                 getMobility(ebosSimulator, perf, mob);
+                const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
+                const double Tw = well_index_[perf] * trans_mult;
                 std::vector<EvalWell> 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<double>& rates) {
+        auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& 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<double>& rates) {
+        auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& 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<typename TypeTag>
+    std::vector<double>
+    MultisegmentWell<TypeTag>::
+    computeCurrentWellRates(const Simulator& ebosSimulator,
+                            DeferredLogger& deferred_logger) const
+    {
+        // Calculate the rates that follow from the current primary variables.
+        std::vector<EvalWell> 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<EvalWell> mob(num_components_, 0.0);
+                getMobility(ebosSimulator, perf, mob);
+                const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
+                const double Tw = well_index_[perf] * trans_mult;
+                std::vector<EvalWell> 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<double> 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

opm/simulators/wells/SimFIBODetails.hpp

@@ -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 <http://www.gnu.org/licenses/>.
-*/
-#ifndef OPM_SIM_FIBO_DETAILS_HPP
-#define OPM_SIM_FIBO_DETAILS_HPP
-
-#include <utility>
-#include <algorithm>
-#include <locale>
-
-#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
-
-namespace Opm
-{
-    namespace SimFIBODetails {
-        typedef std::unordered_map<std::string, Well > 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<double>&            rates)
-        {
-            assert (! p.prediction_mode);
-            assert (rates.size() ==
-                    std::vector<double>::size_type(pu.num_phases));
-
-            if (pu.phase_used[ BlackoilPhases::Aqua ]) {
-                const std::vector<double>::size_type
-                    i = pu.phase_pos[ BlackoilPhases::Aqua ];
-
-                rates[i] = p.water_rate;
-            }
-
-            if (pu.phase_used[ BlackoilPhases::Liquid ]) {
-                const std::vector<double>::size_type
-                    i = pu.phase_pos[ BlackoilPhases::Liquid ];
-
-                rates[i] = p.oil_rate;
-            }
-
-            if (pu.phase_used[ BlackoilPhases::Vapour ]) {
-                const std::vector<double>::size_type
-                    i = pu.phase_pos[ BlackoilPhases::Vapour ];
-
-                rates[i] = p.gas_rate;
-            }
-        }
-    } // namespace SimFIBODetails
-} // namespace Opm
-
-#endif // OPM_SIM_FIBO_DETAILS_HPP

opm/simulators/wells/StandardWell.hpp

@@ -296,6 +296,9 @@ namespace Opm
         using Base::phaseUsage;
         using Base::vfp_properties_;
 
+        virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
+                                                            DeferredLogger& deferred_logger) const override;
+
     protected:
 
         // protected functions from the Base class

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<double>
     StandardWell<TypeTag>::
     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<double>& rates) {
+        auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<double>& 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<double>& rates) {
+        auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& 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<typename TypeTag>
+    std::vector<double>
+    StandardWell<TypeTag>::
+    computeCurrentWellRates(const Simulator& ebosSimulator,
+                            DeferredLogger& deferred_logger) const
+    {
+        // Calculate the rates that follow from the current primary variables.
+        std::vector<EvalWell> 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<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.});
+            getMobility(ebosSimulator, perf, mob, deferred_logger);
+            std::vector<EvalWell> 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<double>(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<double> 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

opm/simulators/wells/WellGroupHelpers.cpp

@@ -22,6 +22,7 @@
 #include <opm/simulators/wells/TargetCalculator.hpp>
 
 #include <algorithm>
+#include <stack>
 #include <vector>
 
 namespace {
@@ -585,6 +586,120 @@ namespace WellGroupHelpers
         wellState.setCurrentInjectionREINRates(group.name(), rein);
     }
 
+
+
+
+    std::map<std::string, double>
+    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<std::string> children;
+        std::set<std::string> leaf_nodes;
+        std::vector<std::string> 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<std::string, std::vector<double>> 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<double>& up = node_inflows[(*upbranch).uptree_node()];
+                const std::vector<double>& 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<std::string, double> 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)
     {

opm/simulators/wells/WellGroupHelpers.hpp

@@ -26,10 +26,13 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/simulators/utils/DeferredLogger.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
+#include <opm/simulators/wells/VFPProdProperties.hpp>
 #include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
 
 #include <algorithm>
 #include <cassert>
+#include <map>
+#include <string>
 #include <type_traits>
 #include <vector>
 
@@ -258,6 +261,11 @@ namespace WellGroupHelpers
                              const WellStateFullyImplicitBlackoil& wellStateNupcol,
                              WellStateFullyImplicitBlackoil& wellState);
 
+    std::map<std::string, double>
+    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);
 

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<double> 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<double> dynamic_thp_limit_;
+
         const PhaseUsage& phaseUsage() const;
 
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const;

opm/simulators/wells/WellInterface_impl.hpp

@@ -328,6 +328,18 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    setDynamicThpLimit(const double thp_limit)
+    {
+       dynamic_thp_limit_ = thp_limit;
+    }
+
+
+
+
+
     template<typename TypeTag>
     double
     WellInterface<TypeTag>::
@@ -403,6 +415,10 @@ namespace Opm
     WellInterface<TypeTag>::
     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<TypeTag>::
     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 <typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    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<double> 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