Merge pull request #2845 from bska/welpi

Add Calculator for Connection/Well Productivity Index

CMakeLists_files.cmake

@@ -44,6 +44,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/VFPProdProperties.cpp
   opm/simulators/wells/VFPInjProperties.cpp
   opm/simulators/wells/WellGroupHelpers.cpp
+  opm/simulators/wells/WellProdIndexCalculator.cpp
   )
 
 if(CUDA_FOUND)
@@ -90,6 +91,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_stoppedwells.cpp
   tests/test_relpermdiagnostics.cpp
   tests/test_norne_pvt.cpp
+  tests/test_wellprodindexcalculator.cpp
   tests/test_wellstatefullyimplicitblackoil.cpp
   )
 
@@ -228,6 +230,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/wells/WellHelpers.hpp
   opm/simulators/wells/WellInterface.hpp
   opm/simulators/wells/WellInterface_impl.hpp
+  opm/simulators/wells/WellProdIndexCalculator.hpp
   opm/simulators/wells/StandardWell.hpp
   opm/simulators/wells/StandardWell_impl.hpp
   opm/simulators/wells/MultisegmentWell.hpp

opm/simulators/wells/WellProdIndexCalculator.cpp

@@ -0,0 +1,131 @@
+/*
+  Copyright 2020 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Connection.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/WellConnections.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+
+namespace {
+    void checkSizeCompatibility(const Opm::WellProdIndexCalculator& wellPICalc,
+                                const std::vector<double>&          connMobility)
+    {
+        if (connMobility.size() != wellPICalc.numConnections()) {
+            throw std::logic_error {
+                "Mobility vector size does not match expected number of connections"
+            };
+        }
+    }
+
+    double logRescale(const double r0, const double rw, const double rd, const double S)
+    {
+        const auto numerator = std::log(r0 / rw) + S;
+        const auto denom     = std::log(rd / rw) + S;
+
+        return numerator / denom;
+    }
+
+    void standardConnFactorsExplicitDrainRadius(const Opm::Well&     well,
+                                                std::vector<double>& stdConnFact)
+    {
+        const auto& connections = well.getConnections();
+        const auto rdrain = well.getDrainageRadius();
+
+        std::transform(connections.begin(), connections.end(), stdConnFact.begin(),
+            [rdrain](const Opm::Connection& conn)
+        {
+            return conn.CF() * logRescale(conn.r0(), conn.rw(), rdrain, conn.skinFactor());
+        });
+    }
+
+    void standardConnFactorsDrainIsEquivalent(const Opm::Well&     well,
+                                              std::vector<double>& stdConnFact)
+    {
+        const auto& connections = well.getConnections();
+
+        std::transform(connections.begin(), connections.end(), stdConnFact.begin(),
+            [](const Opm::Connection& conn)
+        {
+            return conn.CF();
+        });
+    }
+
+    std::vector<double> calculateStandardConnFactors(const Opm::Well& well)
+    {
+        std::vector<double> stdConnFact(well.getConnections().size());
+
+        if (well.getDrainageRadius() > 0.0) {
+            // Well has an explicit drainage radius.  Apply logarithmic
+            // scaling to the CTFs.
+            standardConnFactorsExplicitDrainRadius(well, stdConnFact);
+        }
+        else {
+            // Unspecified drainage radius.  Standard mobility connection
+            // scaling factors are just the regular CTFs.
+            standardConnFactorsDrainIsEquivalent(well, stdConnFact);
+        }
+
+        return stdConnFact;
+    }
+} // namespace Anonymous
+
+Opm::WellProdIndexCalculator::WellProdIndexCalculator(const Well& well)
+    : standardConnFactors_{ calculateStandardConnFactors(well) }
+{}
+
+double
+Opm::WellProdIndexCalculator::
+connectionProdIndStandard(const std::size_t connIdx,
+                          const double      connMobility) const
+{
+    return this->standardConnFactors_[connIdx] * connMobility;
+}
+
+// ===========================================================================
+
+std::vector<double>
+Opm::connectionProdIndStandard(const WellProdIndexCalculator& wellPICalc,
+                               const std::vector<double>&     connMobility)
+{
+    checkSizeCompatibility(wellPICalc, connMobility);
+
+    const auto nConn = wellPICalc.numConnections();
+    auto connPI = connMobility;
+    for (auto connIx = 0*nConn; connIx < nConn; ++connIx) {
+        connPI[connIx] = wellPICalc
+            .connectionProdIndStandard(connIx, connMobility[connIx]);
+    }
+
+    return connPI;
+}
+
+double Opm::wellProdIndStandard(const WellProdIndexCalculator& wellPICalc,
+                                const std::vector<double>&     connMobility)
+{
+    const auto connPI = connectionProdIndStandard(wellPICalc, connMobility);
+
+    return std::accumulate(connPI.begin(), connPI.end(), 0.0);
+}

opm/simulators/wells/WellProdIndexCalculator.hpp

@@ -0,0 +1,105 @@
+/*
+  Copyright 2020 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_WELLPRODINDEXCALCULATOR_HEADER_INCLUDED
+#define OPM_WELLPRODINDEXCALCULATOR_HEADER_INCLUDED
+
+#include <cstddef>
+#include <vector>
+
+namespace Opm {
+    class Well;
+} // namespace Opm
+
+namespace Opm {
+
+    /// Collect per-connection static information to enable calculating
+    /// connection-level or well-level productivity index values when
+    /// incorporating dynamic phase mobilities.
+    class WellProdIndexCalculator
+    {
+    public:
+        /// Constructor
+        ///
+        /// \param[in] well Individual well for which to collect
+        ///   per-connection static data.
+        explicit WellProdIndexCalculator(const Well& well);
+
+        /// Compute connection-level steady-state productivity index value
+        /// using dynamic phase mobility.
+        ///
+        /// \param[in] connIdx Linear connection index.  Must be in the
+        ///   range 0..numConnections() - 1.
+        ///
+        /// \param[in] connMobility Phase mobility at connection \p connIdx.
+        ///   Typically derived from dynamic flow state conditions in cell
+        ///   intersected by well's connection \p connIdx.
+        ///
+        /// \return Connection-level steady-state productivity index.
+        double connectionProdIndStandard(const std::size_t connIdx,
+                                         const double      connMobility) const;
+
+        /// Number of connections in this well.
+        ///
+        /// Used primarily for consistency checks.
+        std::size_t numConnections() const
+        {
+            return this->standardConnFactors_.size();
+        }
+
+    private:
+        /// Static, per-connection multiplicative PI factors.
+        ///
+        /// Corresponds to the well's connection transmissibility factors,
+        /// multiplied by a ratio of logarithms if the well has an explicit,
+        /// positive drainage radius.
+        std::vector<double> standardConnFactors_{};
+    };
+
+    /// Compute connection-level productivity index values for all
+    /// connections in a well.
+    ///
+    /// \param[in] wellPICalc Productivity index calculator.
+    ///
+    /// \param[in] connMobility Phase mobility for each connection.
+    ///   Typically derived from dynamic flow state conditions in cells
+    ///   intersected by well's connections.  Must have one value for each
+    ///   \code wellPICalc.numConnections() \endcode well connection.
+    ///
+    /// \return Connection-level steady-state productivity index values for
+    ///   all connections.
+    std::vector<double>
+    connectionProdIndStandard(const WellProdIndexCalculator& wellPICalc,
+                              const std::vector<double>&     connMobility);
+
+    /// Compute well-level productivity index value.
+    ///
+    /// \param[in] wellPICalc Productivity index calculator.
+    ///
+    /// \param[in] connMobility Phase mobility for each connection.
+    ///   Typically derived from dynamic flow state conditions in cells
+    ///   intersected by well's connections.  Must have one value for each
+    ///   \code wellPICalc.numConnections() \endcode well connection.
+    ///
+    /// \return Well-level steady-state productivity index value.
+    double wellProdIndStandard(const WellProdIndexCalculator& wellPICalc,
+                               const std::vector<double>&     connMobility);
+} // namespace Opm
+
+#endif // OPM_WELLPRODINDEXCALCULATOR_HEADER_INCLUDED

tests/test_wellprodindexcalculator.cpp

@@ -0,0 +1,570 @@
+/*
+  Copyright 2020 Equinor.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE TestWellProdIndexCalculator
+
+#include <boost/test/unit_test.hpp>
+
+#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+#include <opm/parser/eclipse/Units/Units.hpp>
+
+#include <cmath>
+#include <cstddef>
+#include <string>
+
+namespace {
+    double cp_rm3_per_db()
+    {
+        return Opm::prefix::centi*Opm::unit::Poise * Opm::unit::cubic(Opm::unit::meter)
+            / (Opm::unit::day * Opm::unit::barsa);
+    }
+
+    std::string drainRadDefaulted()
+    {
+        return { R"(
+  'P' 'G' 10 10 2005 'LIQ' /
+)" };
+    }
+
+    std::string explicitDrainRad()
+    {
+        // rd = exp(2)
+        return { R"(
+  'P' 'G' 10 10 2005 'LIQ' 7.38905609893065 /
+)" };
+    }
+
+    std::string noSkinFactor_SameCF()
+    {
+        // r0 = exp(1)
+        return { R"(
+  'P' 0 0 1 3 OPEN 1 100 2.0 4* 2.718281828459045 /
+)" };
+    }
+
+    std::string noSkinFactor_DifferentCF()
+    {
+        // r0 = exp(1)
+        return { R"(
+  'P' 0 0 1 1 OPEN 1  50 2.0 4* 2.718281828459045 /
+  'P' 0 0 2 2 OPEN 1 100 2.0 4* 2.718281828459045 /
+  'P' 0 0 3 3 OPEN 1 200 2.0 4* 2.718281828459045 /
+)" };
+    }
+
+    std::string skin2_SameCF()
+    {
+        // r0 = exp(1), Skin = 2
+        return { R"(
+  'P' 0 0 1 3 OPEN 1 100 2.0 1* 2.0 2* 2.718281828459045 /
+)" };
+    }
+
+    std::string skin421_DifferentCF()
+    {
+        // r0 = exp(1), Skin = 4, 2, 1
+        return { R"(
+  'P' 0 0 1 1 OPEN 1  50 2.0 1* 4.0 2* 2.718281828459045 /
+  'P' 0 0 2 2 OPEN 1 100 2.0 1* 2.0 2* 2.718281828459045 /
+  'P' 0 0 3 3 OPEN 1 200 2.0 1* 1.0 2* 2.718281828459045 /
+)" };
+    }
+
+    Opm::Well createWell(const std::string& welspecs,
+                         const std::string& compdat)
+    {
+        const auto deck = Opm::Parser{}.parseString(R"(RUNSPEC
+DIMENS
+  10 10 3 /
+START
+ 8 OCT 2020 /
+GRID
+DXV
+  10*100.0 /
+DYV
+  10*100.0 /
+DZV
+  3*10.0 /
+DEPTHZ
+  121*2000.0 /
+PERMX
+  300*100.0 /
+PERMY
+  300*100.0 /
+PERMZ
+  300*10.0 /
+PORO
+  300*0.3 /
+SCHEDULE
+WELSPECS
+)" + welspecs + R"(
+/
+
+COMPDAT
+)" + compdat + R"(
+/
+
+TSTEP
+  10
+/
+END
+)");
+
+        const auto es    = Opm::EclipseState{ deck };
+        const auto sched = Opm::Schedule{ deck, es };
+
+        return sched.getWell("P", 0);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE(ConnectionLevel)
+
+BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), noSkinFactor_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 1.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 2.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 4.0 * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), noSkinFactor_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 2.0), expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 1.0), expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 0.5), expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), skin2_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 1.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 2.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 4.0 * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), skin421_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 2.0), 1.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 1.0), 1.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 0.5), 1.0 * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), noSkinFactor_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 0.5 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 2.0 * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), noSkinFactor_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 0.25 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.0  * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 4.0  * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), skin2_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 0.75 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.5  * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 3.0  * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), skin421_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), (5.0 / 6.0) * 0.5 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.5         * 1.0 * expectCF, 1.0e-10);
+    BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), (8.0 / 3.0) * 2.0 * expectCF, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_SUITE_END() // ConnectionLevel
+
+// ===========================================================================
+
+BOOST_AUTO_TEST_SUITE(AllConnections)
+
+BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), noSkinFactor_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        1.0*expectCF, 2.0*expectCF, 4.0*expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), noSkinFactor_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 2.0, 1.0, 0.5 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        expectCF, expectCF, expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), skin2_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        1.0*expectCF, 2.0*expectCF, 4.0*expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), skin421_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 2.0, 1.0, 0.5 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        expectCF, expectCF, expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), noSkinFactor_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        0.5*expectCF, 1.0*expectCF, 2.0*expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), noSkinFactor_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        0.25*expectCF, 1.0*expectCF, 4.0*expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), skin2_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        0.75*expectCF, 1.5*expectCF, 3.0*expectCF
+    };
+
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), skin421_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = std::vector<double> {
+        (5.0 / 6.0) * 0.5 * expectCF,
+        1.5         * 1.0 * expectCF,
+        (8.0 / 3.0) * 2.0 * expectCF,
+    };
+    
+    const auto connPI = connectionProdIndStandard(wpiCalc, connMobility);
+    BOOST_CHECK_CLOSE(connPI[0], expectPI[0], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[1], expectPI[1], 1.0e-10);
+    BOOST_CHECK_CLOSE(connPI[2], expectPI[2], 1.0e-10);
+}
+
+BOOST_AUTO_TEST_SUITE_END() // AllConnections
+
+// ===========================================================================
+
+BOOST_AUTO_TEST_SUITE(WellLevel)
+
+BOOST_AUTO_TEST_CASE(allDefaulted_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), noSkinFactor_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (1.0 + 2.0 + 4.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), noSkinFactor_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 2.0, 1.0, 0.5 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (1.0 + 1.0 + 1.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), skin2_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (1.0 + 2.0 + 4.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(drainRadDefaulted(), skin421_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 2.0, 1.0, 0.5 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (1.0 + 1.0 + 1.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), noSkinFactor_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (0.5 + 1.0 + 2.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), noSkinFactor_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (0.25 + 1.0 + 4.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), skin2_SameCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = (0.75 + 1.5 + 3.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF)
+{
+    const auto wpiCalc = Opm::WellProdIndexCalculator {
+        createWell(explicitDrainRad(), skin421_DifferentCF())
+    };
+
+    BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3});
+
+    const auto connMobility = std::vector<double> { 1.0, 2.0, 4.0 };
+
+    const auto expectCF = 100*cp_rm3_per_db();
+    const auto expectPI = ((5.0 / 6.0) * 0.5 +
+                           1.5         * 1.0 +
+                           (8.0 / 3.0) * 2.0)*expectCF;
+
+    BOOST_CHECK_CLOSE(wellProdIndStandard(wpiCalc, connMobility), expectPI, 1.0e-10);
+}
+
+BOOST_AUTO_TEST_SUITE_END() // WellLevel