Use gas lift optimization to get ALQ for well test

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -666,8 +666,11 @@ namespace Opm {
             }
 
             try {
+                GLiftEclWells ecl_well_map;
+                initGliftEclWellMap(ecl_well_map);
                 well->wellTesting(simulator_, simulationTime, this->wellState(),
-                                  this->groupState(), this->wellTestState(), deferred_logger);
+                                  this->groupState(), this->wellTestState(), this->phase_usage_, 
+                                  ecl_well_map, deferred_logger);
             } catch (const std::exception& e) {
                 const std::string msg = fmt::format("Exception during testing of well: {}. The well will not open.\n Exception message: {}", wellEcl.name(), e.what());
                 deferred_logger.warning("WELL_TESTING_FAILED", msg);

opm/simulators/wells/GasLiftSingleWellGeneric.cpp

@@ -170,6 +170,74 @@ runOptimize(const int iteration_idx)
     return state;
 }
 
+template<class Scalar>
+std::pair<Scalar, bool>
+GasLiftSingleWellGeneric<Scalar>::
+wellTestALQ()
+{
+    // If WLIFTOPT item 2 is NO we don't optimize
+    if (!this->optimize_) {
+        return {0.0, false};
+    }
+
+    Scalar temp_alq = std::max(this->min_alq_, 0.0);
+    auto cur_alq = temp_alq;
+    auto init_rates = computeLimitedWellRatesWithALQ_(temp_alq);
+    LimitedRates new_rates = *init_rates;
+    std::optional<LimitedRates> rates;
+    Scalar old_gradient = 0.0;
+    bool alq_is_limited = false;
+    bool increase = true;
+    OptimizeState state {*this, increase};
+    bool success = false;
+
+    while (!state.stop_iteration && (++state.it <= this->max_iterations_)) {
+        if (state.checkAlqOutsideLimits(temp_alq, new_rates.oil))
+            break;
+        std::optional<Scalar> alq_opt;
+        std::tie(alq_opt, alq_is_limited) = state.addOrSubtractAlqIncrement(temp_alq);
+        if (!alq_opt)
+            break;
+
+        temp_alq = *alq_opt;
+        if (this->debug)
+            state.debugShowIterationInfo(temp_alq);
+        rates = new_rates;
+        auto temp_rates = computeLimitedWellRatesWithALQ_(temp_alq);
+        if (!temp_rates)
+            temp_rates->oil = 0.0;
+        if (temp_rates->bhp_is_limited)
+            state.stop_iteration = true;
+
+        auto gradient = state.calcEcoGradient(rates->oil, temp_rates->oil, rates->gas, temp_rates->gas);
+        if (this->debug)
+            debugCheckNegativeGradient_(gradient,
+                                        cur_alq,
+                                        temp_alq,
+                                        rates->oil,
+                                        temp_rates->oil,
+                                        rates->gas,
+                                        temp_rates->gas,
+                                        increase);
+        if (!success && gradient != old_gradient) {
+            success = true;
+        }
+        if (success && gradient <= old_gradient) {
+            break;
+        }
+        cur_alq = temp_alq;
+        new_rates = *temp_rates;
+        old_gradient = gradient;
+    }
+    if (state.it > this->max_iterations_) {
+        warnMaxIterationsExceeded_();
+    }
+    if (success) {
+        this->well_state_.gliftUpdateAlqIncreaseCount(this->well_name_, increase);
+    }
+    return {cur_alq, success};
+}
+
 /****************************************
  * Protected methods in alphabetical order
  ****************************************/

opm/simulators/wells/GasLiftSingleWellGeneric.hpp

@@ -108,6 +108,8 @@ public:
 
     std::unique_ptr<GasLiftWellState<Scalar>> runOptimize(const int iteration_idx);
 
+    std::pair<Scalar, bool> wellTestALQ();
+
     virtual const WellInterfaceGeneric<Scalar>& getWell() const = 0;
 
 protected:

opm/simulators/wells/WellInterface.hpp

@@ -84,7 +84,8 @@ public:
     using GasLiftSingleWell = ::Opm::GasLiftSingleWell<TypeTag>;
     using GLiftOptWells = typename BlackoilWellModel<TypeTag>::GLiftOptWells;
     using GLiftProdWells = typename BlackoilWellModel<TypeTag>::GLiftProdWells;
-    using GLiftWellStateMap =
+    using GLiftEclWells = typename GasLiftGroupInfo<Scalar>::GLiftEclWells;
+    using GLiftWellStateMap = 
         typename BlackoilWellModel<TypeTag>::GLiftWellStateMap;
     using GLiftSyncGroups = typename GasLiftSingleWellGeneric<Scalar>::GLiftSyncGroups;
 
@@ -298,6 +299,8 @@ public:
                      /* const */ WellState<Scalar>& well_state,
                      const GroupState<Scalar>& group_state,
                      WellTestState& welltest_state,
+                     const PhaseUsage& phase_usage,
+                     GLiftEclWells& ecl_well_map,
                      DeferredLogger& deferred_logger);
 
     void checkWellOperability(const Simulator& simulator,
@@ -312,6 +315,9 @@ public:
 
     void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator,
                                              WellState<Scalar>& well_state,
+                                             const GroupState<Scalar>& group_state,
+                                             const PhaseUsage& phase_usage,
+                                             GLiftEclWells& ecl_well_map,
                                              DeferredLogger& deferred_logger);
 
     // check whether the well is operable under the current reservoir condition
@@ -474,6 +480,15 @@ protected:
                              WellState<Scalar>& well_state,
                              const GroupState<Scalar>& group_state,
                              DeferredLogger& deferred_logger);
+    
+
+    template<class GasLiftSingleWell>
+    std::unique_ptr<GasLiftSingleWell> initializeGliftWellTest_(const Simulator& simulator,
+                                                                WellState<Scalar>& well_state,
+                                                                const GroupState<Scalar>& group_state,
+                                                                const PhaseUsage& phase_usage,
+                                                                GLiftEclWells& ecl_well_map,
+                                                                DeferredLogger& deferred_logger);
 
     Eval getPerfCellPressure(const FluidState& fs) const;
 

opm/simulators/wells/WellInterface_impl.hpp

@@ -380,6 +380,8 @@ namespace Opm
                 /* const */ WellState<Scalar>& well_state,
                 const GroupState<Scalar>& group_state,
                 WellTestState& well_test_state,
+                const PhaseUsage& phase_usage,
+                GLiftEclWells& ecl_well_map,
                 DeferredLogger& deferred_logger)
     {
         deferred_logger.info(" well " + this->name() + " is being tested");
@@ -397,7 +399,12 @@ namespace Opm
             const auto report_step = simulator.episodeIndex();
             const auto& glo = schedule.glo(report_step);
             if (glo.active()) {
-                gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
+                gliftBeginTimeStepWellTestUpdateALQ(simulator, 
+                                                    well_state_copy, 
+                                                    group_state,
+                                                    phase_usage, 
+                                                    ecl_well_map, 
+                                                    deferred_logger);
             }
         }
 
@@ -702,13 +709,11 @@ namespace Opm
         bool converged;
         if (has_thp_limit) {
             well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::THP;
-            converged = gliftBeginTimeStepWellTestIterateWellEquations(
-                simulator, dt, well_state, group_state, deferred_logger);
         }
         else {
             well_state.well(this->indexOfWell()).production_cmode = Well::ProducerCMode::BHP;
-            converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
         }
+        converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
         if (converged) {
             deferred_logger.debug("WellTest: Well equation for well " + this->name() +  " converged");
             return true;
@@ -1001,8 +1006,11 @@ namespace Opm
     void
     WellInterface<TypeTag>::
     gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator,
-                          WellState<Scalar>& well_state,
-                          DeferredLogger& deferred_logger)
+                                        WellState<Scalar>& well_state,
+                                        const GroupState<Scalar>& group_state,
+                                        const PhaseUsage& phase_usage,
+                                        GLiftEclWells& ecl_well_map,
+                                        DeferredLogger& deferred_logger)
     {
         const auto& summary_state = simulator.vanguard().summaryState();
         const auto& well_name = this->name();
@@ -1016,28 +1024,42 @@ namespace Opm
         const auto& glo = schedule.glo(report_step_idx);
         if (!glo.has_well(well_name)) {
             const std::string msg = fmt::format(
-                "GLIFT WTEST: Well {} : Gas Lift not activated: "
+                "GLIFT WTEST: Well {} : Gas lift not activated: "
                 "WLIFTOPT is probably missing. Skipping.", well_name);
             deferred_logger.info(msg);
             return;
         }
         const auto& gl_well = glo.well(well_name);
-        auto& max_alq_optional = gl_well.max_rate();
-        Scalar max_alq;
-        if (max_alq_optional) {
-            max_alq = *max_alq_optional;
+
+        // Use gas lift optimization to get ALQ for well test
+        std::unique_ptr<GasLiftSingleWell> glift =
+            initializeGliftWellTest_<GasLiftSingleWell>(simulator,
+                                                        well_state,
+                                                        group_state,
+                                                        phase_usage,
+                                                        ecl_well_map,
+                                                        deferred_logger);
+        auto [max_alq, success] = glift->wellTestALQ();
+        std::string msg;
+        if (success) {
+            well_state.setALQ(well_name, max_alq);
+            msg = fmt::format(
+                "GLIFT WTEST: Well {} : Setting ALQ to optimized value = {}",
+                well_name, max_alq);
         }
         else {
-            const auto& well_ecl = this->wellEcl();
-            const auto& controls = well_ecl.productionControls(summary_state);
-            const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
-            const auto& alq_values = table.getALQAxis();
-            max_alq = alq_values.back();
+            if (!gl_well.use_glo()) {
+                msg = fmt::format(
+                    "GLIFT WTEST: Well {} : Setting ALQ to WLIFTOPT item 3 = {}",
+                    well_name, well_state.getALQ(well_name));
+                
+            }
+            else {
+                msg = fmt::format(
+                    "GLIFT WTEST: Well {} : Gas lift optimization failed, no ALQ set.",
+                    well_name);
+            }
         }
-        well_state.setALQ(well_name, max_alq);
-        const std::string msg = fmt::format(
-            "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
-            well_name, max_alq);
         deferred_logger.info(msg);
     }
 
@@ -1981,4 +2003,48 @@ namespace Opm
         connII[phase_pos] = connIICalc(mt * fs.invB(this->flowPhaseToModelPhaseIdx(phase_pos)).value());
     }
 
+    template<typename TypeTag>
+    template<class GasLiftSingleWell>
+    std::unique_ptr<GasLiftSingleWell> 
+    WellInterface<TypeTag>::
+    initializeGliftWellTest_(const Simulator& simulator,
+                             WellState<Scalar>& well_state,
+                             const GroupState<Scalar>& group_state,
+                             const PhaseUsage& phase_usage,
+                             GLiftEclWells& ecl_well_map,
+                             DeferredLogger& deferred_logger)
+    {
+        // Instantiate group info object (without initialization) since it is needed in GasLiftSingleWell
+        auto& comm = simulator.vanguard().grid().comm();
+        ecl_well_map.try_emplace(this->name(),  &(this->wellEcl()), this->indexOfWell());
+        GasLiftGroupInfo<Scalar> group_info {
+                ecl_well_map,
+                simulator.vanguard().schedule(),
+                simulator.vanguard().summaryState(),
+                simulator.episodeIndex(),
+                simulator.model().newtonMethod().numIterations(),
+                phase_usage,
+                deferred_logger,
+                well_state,
+                group_state,
+                comm,
+                false
+        };
+
+        // Return GasLiftSingleWell object to use the wellTestALQ() function
+        std::set<int> sync_groups;
+        const auto& summary_state = simulator.vanguard().summaryState();
+        return std::make_unique<GasLiftSingleWell>(*this, 
+                                                    simulator, 
+                                                    summary_state,
+                                                    deferred_logger, 
+                                                    well_state, 
+                                                    group_state,
+                                                    group_info, 
+                                                    sync_groups, 
+                                                    comm, 
+                                                    false);
+        
+    }
+
 } // namespace Opm