Merge pull request #2791 from joakim-hove/alq-output

Use gas lift rate value for WALQ output when gas lift rate optimization is in action

opm/parser/eclipse/EclipseState/Schedule/GasLiftOpt.hpp

@@ -267,6 +267,7 @@ public:
     bool active() const;
     bool has_well(const std::string& well) const;
     bool has_group(const std::string& group) const;
+    std::size_t num_wells() const;
 
     static GasLiftOpt serializeObject();
     bool operator==(const GasLiftOpt& other) const;

src/opm/output/eclipse/Summary.cpp

@@ -545,46 +545,69 @@ double efac( const std::vector<std::pair<std::string,double>>& eff_factors, cons
     return (it != eff_factors.end()) ? it->second : 1.0;
 }
 
+inline bool
+has_vfp_table(const Opm::ScheduleState&            sched_state,
+              int vfp_table_number)
+{
+    return sched_state.vfpprod.has(vfp_table_number);
+}
+
+inline Opm::VFPProdTable::ALQ_TYPE
+alq_type(const Opm::ScheduleState&            sched_state,
+         int vfp_table_number)
+{
+    return sched_state.vfpprod(vfp_table_number).getALQType();
+}
+
 inline quantity artificial_lift_quantity( const fn_args& args ) {
     // Note: This function is intentionally supported only at the well level
     // (meaning there's no loop over args.schedule_wells by intention).  Its
     // purpose is to calculate WALQ only.
-    auto alq = quantity { 0.0, measure::identity };
 
+    // Note: in order to determine the correct dimension to use the Summary code
+    // calls the various evaluator functions with a default constructed fn_args
+    // instance. In the case of the WALQ function this does not really work,
+    // because the correct output dimension depends on exactly what physical
+    // quantity is represented by the ALQ - and that again requires quite some
+    // context to determine correctly. The current hack is that if WLIFTOPT is
+    // configured for at least one well we use dimension
+    // measure::gas_surface_rate - otherwise we use measure::identity.
+
+    auto dimension = measure::identity;
+    const auto& glo = args.schedule[args.sim_step].glo();
+    if (glo.num_wells() != 0)
+        dimension = measure::gas_surface_rate;
+
+    auto zero = quantity{0, dimension};
     if (args.schedule_wells.empty()) {
-        return alq;
+        return zero;
     }
 
     const auto* well = args.schedule_wells.front();
     if (well->isInjector()) {
-        return alq;
+        return zero;
     }
 
     auto xwPos = args.wells.find(well->name());
     if ((xwPos == args.wells.end()) ||
         (xwPos->second.dynamicStatus == Opm::Well::Status::SHUT))
     {
-        return alq;
+        return zero;
     }
 
-    alq.value = well->productionControls(args.st).alq_value;
+    const auto& production = well->productionControls(args.st);
+    if (!glo.has_well(well->name()))
+        return { production.alq_value, dimension};
 
-    return alq;
+    const auto& sched_state = args.schedule[args.sim_step];
+    if (alq_type(sched_state, production.vfp_table_number) != Opm::VFPProdTable::ALQ_TYPE::ALQ_GRAT)
+        return zero;
+
+    const double eff_fac = efac(args.eff_factors, well->name());
+    auto alq_rate = eff_fac * xwPos->second.rates.get(rt::alq, production.alq_value);
+    return { alq_rate, dimension };
 }
 
-inline bool
-has_alq_type(const Opm::ScheduleState&            sched_state,
-             const Opm::Well::ProductionControls& pc)
-{
-    return sched_state.vfpprod.has(pc.vfp_table_number);
-}
-
-inline Opm::VFPProdTable::ALQ_TYPE
-alq_type(const Opm::ScheduleState&            sched_state,
-         const Opm::Well::ProductionControls& pc)
-{
-    return sched_state.vfpprod(pc.vfp_table_number).getALQType();
-}
 
 inline quantity glir( const fn_args& args ) {
     if (args.schedule_wells.empty()) {
@@ -607,17 +630,23 @@ inline quantity glir( const fn_args& args ) {
         }
 
         const auto& production = well->productionControls(args.st);
-        if (! has_alq_type(sched_state, production)) {
+        if (! has_vfp_table(sched_state, production.vfp_table_number)) {
             continue;
         }
 
-        const auto thisAlqType = alq_type(sched_state, production);
-        if (thisAlqType != Opm::VFPProdTable::ALQ_TYPE::ALQ_GRAT) {
-            continue;
+        const auto thisAlqType = alq_type(sched_state, production.vfp_table_number);
+        if (thisAlqType == Opm::VFPProdTable::ALQ_TYPE::ALQ_GRAT) {
+            const double eff_fac = efac(args.eff_factors, well->name());
+            alq_rate += eff_fac * xwPos->second.rates.get(rt::alq, production.alq_value);
         }
 
-        const double eff_fac = efac(args.eff_factors, well->name());
-        alq_rate += eff_fac * xwPos->second.rates.get(rt::alq, production.alq_value);
+        if (thisAlqType == Opm::VFPProdTable::ALQ_TYPE::ALQ_IGLR) {
+            const double eff_fac = efac(args.eff_factors, well->name());
+            auto glr = production.alq_value;
+            auto wpr = xwPos->second.rates.get(rt::wat);
+            auto opr = xwPos->second.rates.get(rt::oil);
+            alq_rate += eff_fac * glr * (wpr + opr);
+        }
     }
 
     return { alq_rate, measure::gas_surface_rate };

src/opm/parser/eclipse/EclipseState/Schedule/GasLiftOpt.cpp

@@ -72,6 +72,10 @@ bool GasLiftOpt::has_well(const std::string& wname) const {
     return (iter != this->m_wells.end());
 }
 
+std::size_t GasLiftOpt::num_wells() const {
+    return this->m_wells.size();
+}
+
 bool GasLiftOpt::has_group(const std::string& gname) const {
     const auto iter = this->m_groups.find(gname);
     return (iter != this->m_groups.end());

tests/test_Summary.cpp

@@ -1208,9 +1208,9 @@ BOOST_AUTO_TEST_CASE(GLIR_and_ALQ)
     BOOST_CHECK_CLOSE(1234.56 + 2345.67 + 3456.78,
                       ecl_sum_get_group_var(resp, 1, "B1", "GGLIR"), 1.0e-5);
 
-    BOOST_CHECK_CLOSE(0.0, ecl_sum_get_well_var(resp, 1, "B-1H", "WALQ"), 1.0e-5);
-    BOOST_CHECK_CLOSE(0.0, ecl_sum_get_well_var(resp, 1, "B-2H", "WALQ"), 1.0e-5);
-    BOOST_CHECK_CLOSE(0.0, ecl_sum_get_well_var(resp, 1, "B-3H", "WALQ"), 1.0e-5);
+    BOOST_CHECK_EQUAL(ecl_sum_get_well_var(resp, 1, "B-1H", "WGLIR"), ecl_sum_get_well_var(resp, 1, "B-1H", "WALQ"));
+    BOOST_CHECK_EQUAL(ecl_sum_get_well_var(resp, 1, "B-2H", "WGLIR"), ecl_sum_get_well_var(resp, 1, "B-2H", "WALQ"));
+    BOOST_CHECK_EQUAL(ecl_sum_get_well_var(resp, 1, "B-3H", "WGLIR"), ecl_sum_get_well_var(resp, 1, "B-3H", "WALQ"));
 }
 
 BOOST_AUTO_TEST_CASE(connection_kewords) {