Implement extended network model.

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/wells/BlackoilWellModel.hpp

@@ -196,8 +196,8 @@ namespace Opm {
             {
                 auto grp_nwrk_values = ::Opm::data::GroupAndNetworkValues{};
 
-                this->assignGroupValues(reportStepIdx, sched,
+                this->assignGroupValues(reportStepIdx, sched, grp_nwrk_values.groupData);
-                                        grp_nwrk_values.groupData);
+                this->assignNodeValues(reportStepIdx, sched, grp_nwrk_values.nodeData);
 
                 return grp_nwrk_values;
             }
@@ -315,6 +315,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 +357,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 +455,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

@@ -329,6 +329,8 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     beginTimeStep() {
 
+        updatePerforationIntensiveQuantities();
+
         Opm::DeferredLogger local_deferredLogger;
 
         well_state_ = previous_well_state_;
@@ -392,6 +394,11 @@ namespace Opm {
             local_deferredLogger.warning("WELL_POTENTIAL_CALCULATION_FAILED", msg);
         }
 
+        // 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 +1226,8 @@ namespace Opm {
 
         updateAndCommunicateGroupData();
 
+        updateNetworkPressures();
+
         std::set<std::string> switched_wells;
         std::set<std::string> switched_groups;
 
@@ -1257,6 +1266,44 @@ 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 (producers only, TODO address injectors).
+        for (auto& well : well_container_) {
+            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);
+                }
+            }
+        }
+
+        // Output debug info.
+        if (terminal_output_) {
+            std::ostringstream oss;
+            oss << "Node pressures in network nodes, in bar:\n";
+            for (const auto& [name, pressure] : node_pressures_) {
+                oss << name << "  " << unit::convert::to(pressure, unit::barsa) << '\n';
+            }
+            OpmLog::debug(oss.str());
+        }
+    }
+
+
+
 
     template<typename TypeTag>
     void
@@ -2503,6 +2550,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_impl.hpp

@@ -2026,13 +2026,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);
@@ -3306,11 +3306,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 +3530,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.

opm/simulators/wells/StandardWell.hpp

@@ -296,6 +296,10 @@ namespace Opm
         using Base::phaseUsage;
         using Base::vfp_properties_;
 
+        virtual void updateWellStateRates(const Simulator& ebosSimulator,
+                                          WellState& well_state,
+                                          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,
+                                    const SummaryState& summary_state,
-                                   DeferredLogger& deferred_logger,
+                                    DeferredLogger& deferred_logger,
-                                   double alq_value) const
+                                    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,85 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    void
+    StandardWell<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<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];
+            }
+        }
+
+
+        // We must keep in mind that component and phase indices are different here.
+        // Therefore we set up a mapping to make the last code block simpler.
+        const auto& pu = phaseUsage();
+        const int wpi = Water;
+        const int opi = Oil;
+        const int gpi = Gas;
+        const unsigned int wci = FluidSystem::waterCompIdx;
+        const unsigned int oci = FluidSystem::oilCompIdx;
+        const unsigned int gci = FluidSystem::gasCompIdx;
+        std::pair<int, unsigned int> phase_comp[3] = { {wpi, wci},
+                                                       {opi, oci},
+                                                       {gpi, gci} };
+        std::vector<int> phase_to_comp(number_of_phases_, -1);
+        for (const auto& pc : phase_comp) {
+            if (pu.phase_used[pc.first]) {
+                const int phase_idx = pu.phase_pos[pc.first];
+                const int comp_idx = Indices::canonicalToActiveComponentIndex(pc.second);
+                phase_to_comp[phase_idx] = comp_idx;
+            }
+        }
+
+        // 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 = phase_to_comp[nonzero_rate_index];
+        for (int p = 0; p < number_of_phases_; ++p) {
+            if (p != nonzero_rate_index) {
+                const int comp_idx = phase_to_comp[p];
+                double& rate = well_state.wellRates()[index_of_well_ * number_of_phases_ + p];
+                rate = (initial_nonzero_rate/well_q_s[comp_idx_nz].value()) * (well_q_s[comp_idx].value());
+            }
+        }
+    }
+
+
 } // 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 1
+#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,14 @@ namespace Opm
         // update perforation water throughput based on solved water rate
         virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
 
+        virtual void updateWellStateRates(const Simulator& ebosSimulator,
+                                          WellState& well_state,
+                                          DeferredLogger& deferred_logger) const
+        {
+            // Default: do nothing.
+            // TODO: make this a pure virtual.
+        }
+
         void stopWell() {
             wellIsStopped_ = true;
         }
@@ -288,6 +296,7 @@ namespace Opm
 
         void setWsolvent(const double wsolvent);
 
+        void setDynamicThpLimit(const double thp_limit);
 
     protected:
 
@@ -384,6 +393,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;