From d2259017333aa6895131caa68a3a64ef8c286248 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Wed, 20 Sep 2023 22:54:48 +0200
Subject: [PATCH] Add Support for WPIx Summary Vectors

We calculate these, for now, as

   WPIx = WOPR / (WBPx - WBHP)

but it is possible that we may have to revise this definition,
especially for injectors.
---
 src/opm/output/eclipse/Summary.cpp | 41 ++++++++++++++++++++++++++++++
 tests/summary_deck.DATA            | 13 ++++++----
 tests/test_Summary.cpp             |  7 +++++
 3 files changed, 56 insertions(+), 5 deletions(-)

diff --git a/src/opm/output/eclipse/Summary.cpp b/src/opm/output/eclipse/Summary.cpp
index d79a4e422..59f60798f 100644
--- a/src/opm/output/eclipse/Summary.cpp
+++ b/src/opm/output/eclipse/Summary.cpp
@@ -1535,6 +1535,41 @@ inline quantity potential_rate( const fn_args& args )
     return { sum, rate_unit< phase >() };
 }
 
+template <Opm::data::WellBlockAvgPress::Quantity wbp_quantity>
+quantity well_block_average_prod_index(const fn_args& args)
+{
+    // Note: This WPIn evaluation function is supported only at the well
+    // level.  There is intentionally no loop over args.schedule_wells.
+
+    const auto unit = rate_unit<rt::productivity_index_oil>();
+    const quantity zero = { 0.0, unit };
+
+    if (args.schedule_wells.empty()) {
+        return zero;
+    }
+
+    const auto& name = args.schedule_wells.front()->name();
+
+    auto xwPos = args.wells.find(name);
+    if ((xwPos == args.wells.end()) ||
+        (xwPos->second.dynamicStatus == Opm::Well::Status::SHUT))
+    {
+        return zero;
+    }
+
+    auto p = args.wbp.values.find(args.schedule_wells.front()->name());
+    if (p == args.wbp.values.end()) {
+        return zero;
+    }
+
+    // Rt::oil is intentional.
+    const auto eff_fac = efac(args.eff_factors, name);
+    const auto q  = xwPos->second.rates.get(rt::oil, 0.0) * eff_fac;
+    const auto dp = p->second[wbp_quantity] - xwPos->second.bhp;
+
+    return { - q / dp, unit };
+}
+
 inline quantity preferred_phase_productivty_index(const fn_args& args)
 {
     if (args.schedule_wells.empty())
@@ -2402,6 +2437,12 @@ static const auto funs = std::unordered_map<std::string, ofun> {
     { "WPIG", potential_rate< rt::productivity_index_gas >},
     { "WPIL", sum( potential_rate< rt::productivity_index_water, true, false >,
                    potential_rate< rt::productivity_index_oil, true, false >)},
+
+    { "WPI1", well_block_average_prod_index<Opm::data::WellBlockAvgPress::Quantity::WBP>  },
+    { "WPI4", well_block_average_prod_index<Opm::data::WellBlockAvgPress::Quantity::WBP4> },
+    { "WPI5", well_block_average_prod_index<Opm::data::WellBlockAvgPress::Quantity::WBP5> },
+    { "WPI9", well_block_average_prod_index<Opm::data::WellBlockAvgPress::Quantity::WBP9> },
+
     // Well potential
     { "WWPP", potential_rate< rt::well_potential_water , true, false>},
     { "WOPP", potential_rate< rt::well_potential_oil , true, false>},
diff --git a/tests/summary_deck.DATA b/tests/summary_deck.DATA
index 559809637..bc0a9fe62 100644
--- a/tests/summary_deck.DATA
+++ b/tests/summary_deck.DATA
@@ -716,11 +716,14 @@ WPIW
 /
 WPIL
 /
-WBP
-/
-WBP4
-/
-
+WPI1
+ 'W_1' 'W_2' /
+WPI4
+ 'W_1' /
+WPI5
+ 'W_1' /
+WPI9
+ 'W_1' /
 -- from model2
 WOPTF
 /
diff --git a/tests/test_Summary.cpp b/tests/test_Summary.cpp
index c92c89739..85905f1cc 100644
--- a/tests/test_Summary.cpp
+++ b/tests/test_Summary.cpp
@@ -663,6 +663,13 @@ BOOST_AUTO_TEST_CASE(well_keywords)
     BOOST_CHECK_CLOSE( 10.11, ecl_sum_get_well_var( resp, 1, "W_1", "WPI"  ), 1.0e-5 );
     BOOST_CHECK_CLOSE( 21.01, ecl_sum_get_well_var( resp, 1, "W_1", "WPIL" ), 1.0e-5 );
 
+    BOOST_CHECK_CLOSE(10.1 / (123.456 - 0.1), ecl_sum_get_well_var( resp, 1, "W_1", "WPI1" ), 1.0e-5 );
+    BOOST_CHECK_CLOSE(10.1 / (123.567 - 0.1), ecl_sum_get_well_var( resp, 1, "W_1", "WPI4" ), 1.0e-5 );
+    BOOST_CHECK_CLOSE(10.1 / (123.678 - 0.1), ecl_sum_get_well_var( resp, 1, "W_1", "WPI5" ), 1.0e-5 );
+    BOOST_CHECK_CLOSE(10.1 / (123.789 - 0.1), ecl_sum_get_well_var( resp, 1, "W_1", "WPI9" ), 1.0e-5 );
+
+    BOOST_CHECK_CLOSE(0.0, ecl_sum_get_well_var( resp, 1, "W_2", "WPI1" ), 1.0e-5 );
+
     BOOST_CHECK_CLOSE( 20.9 , ecl_sum_get_well_var( resp, 1, "W_2", "WPIW" ), 1.0e-5 );
     BOOST_CHECK_CLOSE( 20.11, ecl_sum_get_well_var( resp, 1, "W_2", "WPIO" ), 1.0e-5 );
     BOOST_CHECK_CLOSE( 20.12, ecl_sum_get_well_var( resp, 1, "W_2", "WPIG" ), 1.0e-5 );