Merge pull request #4935 from plgbrts/auto-choke

WIP: Enabling automatic choke to model a sub-sea manifold in an extended network

opm/simulators/wells/BlackoilWellModel.hpp

@@ -413,7 +413,7 @@ template<class Scalar> class WellContributions;
             Scalar gravity_{};
             std::vector<Scalar> depth_{};
             bool alternative_well_rate_init_{};
-
+            std::map<std::string, Scalar> well_group_thp_calc_;
             std::unique_ptr<RateConverterType> rateConverter_{};
             std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};
 
@@ -465,6 +465,8 @@ template<class Scalar> class WellContributions;
                                               const double dt,
                                               DeferredLogger& local_deferredLogger);
 
+            void computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger);
+
             /// Update rank's notion of intersecting wells and their
             /// associate solution variables.
             ///

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -28,6 +28,7 @@
 #endif
 
 #include <opm/grid/utility/cartesianToCompressed.hpp>
+#include <opm/common/utility/numeric/RootFinders.hpp>
 
 #include <opm/input/eclipse/Schedule/Network/Balance.hpp>
 #include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
@@ -40,6 +41,9 @@
 #include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
 #include <opm/simulators/wells/ParallelWBPCalculation.hpp>
 #include <opm/simulators/wells/VFPProperties.hpp>
+#include <opm/simulators/wells/WellBhpThpCalculator.hpp>
+#include <opm/simulators/wells/WellGroupHelpers.hpp>
+#include <opm/simulators/wells/TargetCalculator.hpp>
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/utils/MPIPacker.hpp>
@@ -1243,7 +1247,142 @@ namespace Opm {
         return {more_network_update, well_group_control_changed};
     }
 
+    // This function is to be used for well groups in an extended network that act as a subsea manifold
+    // The wells of such group should have a common THP and total phase rate(s) obeying (if possible) 
+    // the well group constraint set by GCONPROD
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger)
+    {
+        const int reportStepIdx = this->simulator_.episodeIndex();
+        const auto& network = this->schedule()[reportStepIdx].network();
+        const auto& balance = this->schedule()[reportStepIdx].network_balance();
+        const Scalar thp_tolerance = balance.thp_tolerance();
 
+        if (!network.active()) {
+            return;
+        }
+
+        auto& well_state = this->wellState();
+        auto& group_state = this->groupState();
+
+        for (const std::string& nodeName : network.node_names()) {
+            const bool has_choke = network.node(nodeName).as_choke();
+            if (has_choke) {
+                const auto& summary_state = this->simulator_.vanguard().summaryState();
+                const Group& group = this->schedule().getGroup(nodeName, reportStepIdx);
+                const auto ctrl = group.productionControls(summary_state);
+                const auto cmode = ctrl.cmode;
+                const auto pu = this->phase_usage_;
+                //TODO: Auto choke combined with RESV control is not supported
+                std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
+                Scalar gratTargetFromSales = 0.0; 
+                if (group_state.has_grat_sales_target(group.name()))
+                    gratTargetFromSales = group_state.grat_sales_target(group.name());
+
+                WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
+                                                  gratTargetFromSales, nodeName, group_state,
+                                                  group.has_gpmaint_control(cmode));
+                const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
+
+                auto mismatch = [&] (auto group_thp) {
+                    Scalar group_rate(0.0);
+                    Scalar rate(0.0);
+                    for (auto& well : this->well_container_) {
+                        std::string well_name = well->name();
+                        auto& ws = well_state.well(well_name);
+                        if (group.hasWell(well_name)) {
+                            well->setDynamicThpLimit(group_thp);
+                            const Well& well_ecl = this->wells_ecl_[well->indexOfWell()];
+                            const auto inj_controls = Well::InjectionControls(0);
+                            const auto prod_controls = well_ecl.productionControls(summary_state);
+                            well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger,  false, false);
+                            rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
+                            group_rate += rate;
+                        }
+                    }
+                    return (group_rate - orig_target)/orig_target;
+                };
+
+                const auto upbranch = network.uptree_branch(nodeName);
+                const auto it = this->node_pressures_.find((*upbranch).uptree_node());
+                const Scalar nodal_pressure = it->second;
+                Scalar well_group_thp = nodal_pressure;
+
+                std::optional<Scalar> autochoke_thp;
+                if (auto iter = this->well_group_thp_calc_.find(nodeName); iter != this->well_group_thp_calc_.end()) {
+                    autochoke_thp = this->well_group_thp_calc_.at(nodeName);
+                }
+
+                //Find an initial bracket
+                std::array<Scalar, 2> range_initial;
+                if (!autochoke_thp.has_value()){
+                    Scalar min_thp, max_thp;
+                    // Retrieve the terminal pressure of the associated root of the manifold group
+                    std::string node_name =  nodeName;
+                    while (!network.node(node_name).terminal_pressure().has_value()) {
+                        auto branch = network.uptree_branch(node_name).value();
+                        node_name = branch.uptree_node();
+                    }
+                    min_thp = network.node(node_name).terminal_pressure().value();
+                    WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp);
+                    // Narrow down the bracket
+                    Scalar low1, high1;
+                    std::array<Scalar, 2> range = {0.9*min_thp, 1.1*max_thp};
+                    std::optional<Scalar> appr_sol;
+                    WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger);
+                    min_thp = low1;
+                    max_thp = high1;
+                    range_initial = {min_thp, max_thp};
+                }
+
+                if (!autochoke_thp.has_value() || autochoke_thp.value() > nodal_pressure) {
+                    // The bracket is based on the initial bracket or on a range based on a previous calculated group thp
+                    std::array<Scalar, 2> range = autochoke_thp.has_value() ?
+                        std::array<Scalar, 2>{0.9 * autochoke_thp.value(), 1.1 * autochoke_thp.value()} : range_initial;
+                    Scalar low, high;
+                    std::optional<Scalar> approximate_solution;
+                    const Scalar tolerance1 = thp_tolerance;
+                    local_deferredLogger.debug("Using brute force search to bracket the group THP");
+                    const bool finding_bracket = WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger);
+
+                    if (approximate_solution.has_value()) {
+                        autochoke_thp = *approximate_solution;
+                        local_deferredLogger.debug("Approximate group THP value found: "  + std::to_string(autochoke_thp.value()));
+                    } else if (finding_bracket) {
+                        const Scalar tolerance2 = thp_tolerance;
+                        const int max_iteration_solve = 100;
+                        int iteration = 0;
+                        autochoke_thp = RegulaFalsiBisection<ThrowOnError>::
+                                         solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration);
+                        local_deferredLogger.debug(" bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " +
+                                                   "iteration = " + std::to_string(iteration));
+                        local_deferredLogger.debug("Group THP value = " + std::to_string(autochoke_thp.value()));
+                    } else {
+                        autochoke_thp.reset();
+                        local_deferredLogger.debug("Group THP solve failed due to bracketing failure");
+                    }
+                }
+                 if (autochoke_thp.has_value()) {
+                    well_group_thp_calc_[nodeName] = autochoke_thp.value();
+                    // Note: The node pressure of the auto-choke node is set to well_group_thp in computeNetworkPressures()
+                    // and must be larger or equal to the pressure of the uptree node of its branch.
+                    well_group_thp = std::max(autochoke_thp.value(), nodal_pressure);
+                }
+
+                for (auto& well : this->well_container_) {
+                    std::string well_name = well->name();
+                    if (group.hasWell(well_name)) {
+                        well->setDynamicThpLimit(well_group_thp);
+                    }
+                }
+
+                // Use the group THP in computeNetworkPressures().
+                group_state.update_well_group_thp(nodeName, well_group_thp);
+            }
+        }
+    }
 
     template<typename TypeTag>
     void
@@ -1905,6 +2044,9 @@ namespace Opm {
         // network related
         bool more_network_update = false;
         if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
+            const double dt = this->simulator_.timeStepSize();
+            // Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES)
+            computeWellGroupThp(dt, deferred_logger);
             const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx);
             const Scalar network_imbalance = comm.max(local_network_imbalance);
             const auto& balance = this->schedule()[episodeIdx].network_balance();

opm/simulators/wells/GroupState.cpp

@@ -102,6 +102,26 @@ GroupState<Scalar>::production_rates(const std::string& gname) const
 
 //-------------------------------------------------------------------------
 
+template<class Scalar>
+void GroupState<Scalar>::
+GroupState::update_well_group_thp(const std::string& gname, const double& thp) 
+{
+    this->group_thp[gname] = thp;
+}
+
+template<class Scalar>
+Scalar GroupState<Scalar>::
+GroupState::well_group_thp(const std::string& gname) const 
+{
+    auto group_iter = this->group_thp.find(gname);
+    if (group_iter == this->group_thp.end())
+        throw std::logic_error("No such group");
+
+    return group_iter->second;
+}
+
+//-------------------------------------------------------------------------
+
 template<class Scalar>
 bool GroupState<Scalar>::
 has_production_reduction_rates(const std::string& gname) const

opm/simulators/wells/GroupState.hpp

@@ -49,6 +49,9 @@ public:
                                  const std::vector<Scalar>& rates);
     const std::vector<Scalar>& production_rates(const std::string& gname) const;
 
+    void update_well_group_thp(const std::string& gname, const double& thp);
+    Scalar well_group_thp(const std::string& gname) const;
+
     bool has_production_reduction_rates(const std::string& gname) const;
     void update_production_reduction_rates(const std::string& gname,
                                            const std::vector<Scalar>& rates);
@@ -175,6 +178,7 @@ public:
         serializer(num_phases);
         serializer(m_production_rates);
         serializer(production_controls);
+        serializer(group_thp);
         serializer(prod_red_rates);
         serializer(inj_red_rates);
         serializer(inj_surface_rates);
@@ -199,6 +203,7 @@ private:
     std::map<std::string, Scalar> inj_vrep_rate;
     std::map<std::string, Scalar> m_grat_sales_target;
     std::map<std::string, Scalar> m_gpmaint_target;
+    std::map<std::string, Scalar> group_thp;
 
     std::map<std::pair<Phase, std::string>, Group::InjectionCMode> injection_controls;
     WellContainer<GPMaint::State> gpmaint_state;

opm/simulators/wells/WellAssemble.cpp

@@ -35,6 +35,7 @@
 #include <opm/simulators/wells/WellHelpers.hpp>
 #include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
 #include <opm/simulators/wells/WellState.hpp>
+#include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
 
 #include <cassert>
 #include <stdexcept>
@@ -162,12 +163,13 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
                 well_.rateConverter().calcCoeff(id, region, coeff);
 
         };
-        WellGroupControls(well_).getGroupProductionControl(group,
+
+        WellGroupControls(well_).getGroupProductionControl(group, 
                                                            well_state,
                                                            group_state,
                                                            schedule,
                                                            summaryState,
-                                                           bhp,
+                                                           bhp, 
                                                            active_rates,
                                                            rCoeff,
                                                            efficiencyFactor,

opm/simulators/wells/WellBhpThpCalculator.cpp

@@ -1005,6 +1005,70 @@ getFloIPR(const WellState<Scalar>& well_state,
                           detail::getFlo(table, aqua_b, liquid_b, vapour_b));
 }
 
+template<class Scalar>
+bool
+WellBhpThpCalculator<Scalar>::
+bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
+                  const std::array<Scalar, 2>& range,
+                  Scalar& low, Scalar& high,
+                  std::optional<Scalar>& approximate_solution,
+                  const Scalar& limit,
+                  DeferredLogger& deferred_logger)
+{
+    bool bracket_found = false;
+    low = range[0];
+    high = range[1];
+    const int sample_number = 300;
+    const Scalar interval = (high - low) / sample_number;
+    Scalar eq_low = eq(low);
+    Scalar eq_high = 0.0;
+    for (int i = 0; i < sample_number + 1; ++i) {
+        high = range[0] + interval * i;
+        eq_high = eq(high);
+        if ( (std::fabs(eq_high) < limit)) {
+            approximate_solution = high;
+            break;
+        }
+        if (eq_high * eq_low <= 0.) {
+            bracket_found = true;
+            break;
+        }
+        low = high;
+        eq_low = eq_high;
+    }
+
+    if (bracket_found) {
+        deferred_logger.debug(
+                " brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) +
+                " high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high));
+    }
+    return bracket_found;
+}
+
+template<class Scalar>
+bool
+WellBhpThpCalculator<Scalar>::
+bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
+                  Scalar& min_thp, Scalar& max_thp)
+{
+    bool bracket_found = false;
+    constexpr int sample_number = 1000; 
+    constexpr Scalar interval = 1E5; 
+    Scalar eq_low = eq(min_thp);
+    Scalar eq_high = 0.0;
+    for (int i = 0; i < sample_number + 1; ++i) {
+        max_thp = min_thp + interval * i;
+        eq_high = eq(max_thp);
+        if (eq_high * eq_low <= 0.) {
+            bracket_found = true;
+            min_thp = max_thp - interval;
+            break;
+        }
+        eq_low = eq_high;
+    }
+    return bracket_found;
+}
+
 template class WellBhpThpCalculator<double>;
 
 #define INSTANCE(...) \

opm/simulators/wells/WellBhpThpCalculator.hpp

@@ -117,7 +117,19 @@ public:
   std::pair<Scalar, Scalar>
   getFloIPR(const WellState<Scalar>& well_state,
             const Well& well, 
-            const SummaryState& summary_state) const;                                                                                              
+            const SummaryState& summary_state) const;
+
+  //! \brief Find limits using brute-force solver.
+  static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
+                                const std::array<Scalar, 2>& range,
+                                Scalar& low, Scalar& high,
+                                std::optional<Scalar>& approximate_solution,
+                                const Scalar& limit,
+                                DeferredLogger& deferred_logger);
+
+  //! \brief Find limits using brute-force solver.
+  static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
+                                Scalar& min_thp, Scalar& max_thp);
 
 private:
     //! \brief Compute BHP from THP limit for an injector - implementation.
@@ -155,11 +167,12 @@ private:
                        std::optional<Scalar>& approximate_solution,
                        DeferredLogger& deferred_logger) const;
 
-    //! \brief Find limits using brute-force solver.
+ //! \brief Find limits using brute-force solver.
-    static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
+  static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
-                                  const std::array<Scalar, 2>& range,
+                                const std::array<Scalar, 2>& range,
-                                  Scalar& low, Scalar& high,
+                                Scalar& low, Scalar& high,
-                                  DeferredLogger& deferred_logger);
+                                DeferredLogger& deferred_logger);
+
 
     Scalar findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
                                      const Scalar bhp,

opm/simulators/wells/WellGroupHelpers.cpp

@@ -41,6 +41,7 @@
 #include <opm/simulators/wells/TargetCalculator.hpp>
 #include <opm/simulators/wells/VFPProdProperties.hpp>
 #include <opm/simulators/wells/WellState.hpp>
+#include <opm/simulators/wells/GroupState.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -935,7 +936,12 @@ computeNetworkPressures(const Network::ExtNetwork& network,
 #endif
                 } else {
                     // Table number specified as 9999 in the deck, no pressure loss.
-                    node_pressures[node] = up_press;
+                    if (network.node(node).as_choke()){
+                        // Node pressure is set to the group THP.
+                        node_pressures[node] = group_state.well_group_thp(node);
+                    } else {
+                            node_pressures[node] = up_press;
+                    }
                 }
             }
         }

opm/simulators/wells/WellInterface.hpp

@@ -374,6 +374,15 @@ public:
     void updateConnectionTransmissibilityFactor(const Simulator& simulator,
                                                 SingleWellState<Scalar>& ws) const;
 
+    virtual bool iterateWellEqWithSwitching(const Simulator& simulator,
+                                            const double dt,
+                                            const WellInjectionControls& inj_controls,
+                                            const WellProductionControls& prod_controls,
+                                            WellState<Scalar>& well_state,
+                                            const GroupState<Scalar>& group_state,
+                                            DeferredLogger& deferred_logger, 
+                                            const bool fixed_control = false, 
+                                            const bool fixed_status = false) = 0;
 protected:
     // simulation parameters
     const ModelParameters& param_;
@@ -428,16 +437,6 @@ protected:
                                           const GroupState<Scalar>& group_state,
                                           DeferredLogger& deferred_logger) = 0;
 
-    virtual bool iterateWellEqWithSwitching(const Simulator& simulator,
-                                            const double dt,
-                                            const WellInjectionControls& inj_controls,
-                                            const WellProductionControls& prod_controls,
-                                            WellState<Scalar>& well_state,
-                                            const GroupState<Scalar>& group_state,
-                                            DeferredLogger& deferred_logger, 
-                                            const bool fixed_control = false, 
-                                            const bool fixed_status = false) = 0;
-
     virtual void updateIPRImplicit(const Simulator& simulator,
                                    WellState<Scalar>& well_state,
                                    DeferredLogger& deferred_logger) = 0;                                            

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -103,6 +103,7 @@ public:
     void setWsolvent(const Scalar wsolvent);
     void setDynamicThpLimit(const Scalar thp_limit);
     std::optional<Scalar> getDynamicThpLimit() const;
+    void setDynamicThpLimit(const std::optional<Scalar> thp_limit);
     void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
 
     /// Returns true if the well has one or more THP limits/constraints.