Merge pull request #2846 from bska/calc-wellpi-endoftimestep

Unconditionally Calculate PI at End of Timestep

opm/simulators/wells/BlackoilWellModel.hpp

@@ -53,6 +53,7 @@
 #include <opm/simulators/wells/StandardWell.hpp>
 #include <opm/simulators/wells/MultisegmentWell.hpp>
 #include <opm/simulators/wells/WellGroupHelpers.hpp>
+#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
@@ -268,6 +269,7 @@ namespace Opm {
 
             std::vector< Well > wells_ecl_;
             std::vector< std::vector<PerforationData> > well_perf_data_;
+            std::vector< WellProdIndexCalculator > prod_index_calc_;
 
             bool wells_active_;
 
@@ -281,6 +283,7 @@ namespace Opm {
 
             std::function<bool(const Well&)> is_shut_or_defunct_;
 
+            void initializeWellProdIndCalculators();
             void initializeWellPerfData();
 
             // create the well container
@@ -377,6 +380,8 @@ namespace Opm {
 
             void calculateEfficiencyFactors(const int reportStepIdx);
 
+            void calculateProductivityIndexValues(DeferredLogger& deferred_logger);
+
             // it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation()
             // makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions
             // twice at the beginning of the time step

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -237,6 +237,8 @@ namespace Opm {
         int globalNumWells = 0;
         // Make wells_ecl_ contain only this partition's non-shut wells.
         wells_ecl_ = getLocalNonshutWells(timeStepIdx, globalNumWells);
+
+        this->initializeWellProdIndCalculators();
         initializeWellPerfData();
 
         // Wells are active if they are active wells on at least
@@ -507,6 +509,8 @@ namespace Opm {
         const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
         checkGconsaleLimits(fieldGroup, well_state_, local_deferredLogger);
 
+        this->calculateProductivityIndexValues(local_deferredLogger);
+
         previous_well_state_ = well_state_;
 
         Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
@@ -566,6 +570,7 @@ namespace Opm {
         // Make wells_ecl_ contain only this partition's non-shut wells.
         wells_ecl_ = getLocalNonshutWells(report_step, globalNumWells);
 
+        this->initializeWellProdIndCalculators();
         initializeWellPerfData();
 
         const int nw = wells_ecl_.size();
@@ -586,6 +591,22 @@ namespace Opm {
 
 
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    initializeWellProdIndCalculators()
+    {
+        this->prod_index_calc_.clear();
+        this->prod_index_calc_.reserve(this->wells_ecl_.size());
+        for (const auto& well : this->wells_ecl_) {
+            this->prod_index_calc_.emplace_back(well);
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -1439,6 +1460,31 @@ namespace Opm {
 
 
 
+
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    calculateProductivityIndexValues(DeferredLogger& deferred_logger)
+    {
+        if (! this->localWellsActive()) {
+            return;
+        }
+
+        auto piCalc = this->prod_index_calc_.begin();
+        for (const auto& wellPtr : this->well_container_) {
+            wellPtr->updateProductivityIndex(this->ebosSimulator_,
+                                             *piCalc,
+                                             this->well_state_,
+                                             deferred_logger);
+            ++piCalc;
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/MultisegmentWell.hpp

@@ -22,9 +22,10 @@
 #ifndef OPM_MULTISEGMENTWELL_HEADER_INCLUDED
 #define OPM_MULTISEGMENTWELL_HEADER_INCLUDED
 
-
 #include <opm/simulators/wells/WellInterface.hpp>
 
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
+
 namespace Opm
 {
 
@@ -43,7 +44,7 @@ namespace Opm
         using typename Base::Indices;
         using typename Base::RateConverterType;
         using typename Base::SparseMatrixAdapter;
-
+        using typename Base::FluidState;
 
         /// the number of reservior equations
         using Base::numEq;
@@ -168,6 +169,11 @@ namespace Opm
                                                  const WellState& well_state,
                                                  Opm::DeferredLogger& deferred_logger) override; // should be const?
 
+        virtual void updateProductivityIndex(const Simulator& ebosSimulator,
+                                             const WellProdIndexCalculator& wellPICalc,
+                                             WellState& well_state,
+                                             DeferredLogger& deferred_logger) const override;
+
         virtual void  addWellContributions(SparseMatrixAdapter& jacobian) const override;
 
         /// number of segments for this well
@@ -179,6 +185,17 @@ namespace Opm
         virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
                                                             DeferredLogger& deferred_logger) const override;
 
+        void computeConnLevelProdInd(const FluidState& fs,
+                                     const std::function<double(const double)>& connPICalc,
+                                     const std::vector<EvalWell>& mobility,
+                                     double* connPI) const;
+
+        void computeConnLevelInjInd(const FluidState& fs,
+                                    const Phase preferred_phase,
+                                    const std::function<double(const double)>& connIICalc,
+                                    const std::vector<EvalWell>& mobility,
+                                    double* connII,
+                                    DeferredLogger& deferred_logger) const;
     protected:
         int number_segments_;
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -23,6 +23,8 @@
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>
 
+#include <algorithm>
+
 namespace Opm
 {
 
@@ -1176,6 +1178,80 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    updateProductivityIndex(const Simulator& ebosSimulator,
+                            const WellProdIndexCalculator& wellPICalc,
+                            WellState& well_state,
+                            DeferredLogger& deferred_logger) const
+    {
+        auto fluidState = [&ebosSimulator, this](const int perf)
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            return ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+        };
+
+        const int np = this->number_of_phases_;
+        auto setToZero = [np](double* x) -> void
+        {
+            std::fill_n(x, np, 0.0);
+        };
+
+        auto addVector = [np](const double* src, double* dest) -> void
+        {
+            std::transform(src, src + np, dest, dest, std::plus<>{});
+        };
+
+        auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
+        auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
+
+        setToZero(wellPI);
+
+        const auto preferred_phase = this->well_ecl_.getPreferredPhase();
+
+        const auto& allConn = this->well_ecl_.getConnections();
+        const auto  nPerf   = allConn.size();
+        auto subsetPerfID   = 0*nPerf;
+        for (auto allPerfID = 0*nPerf; allPerfID < nPerf; ++allPerfID) {
+            if (allConn[allPerfID].state() == Connection::State::SHUT) {
+                continue;
+            }
+
+            auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
+            {
+                return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
+            };
+
+            std::vector<EvalWell> mob(this->num_components_, 0.0);
+            getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob);
+
+            const auto& fs = fluidState(subsetPerfID);
+            setToZero(connPI);
+
+            if (this->isInjector()) {
+                this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
+                                             mob, connPI, deferred_logger);
+            }
+            else {  // Production or zero flow rate
+                this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
+            }
+
+            addVector(connPI, wellPI);
+
+            ++subsetPerfID;
+            connPI += np;
+        }
+
+        assert (static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
+                "Internal logic error in processing connections for PI/II");
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
@@ -3893,4 +3969,80 @@ namespace Opm
 
 
 
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    computeConnLevelProdInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
+                            const std::function<double(const double)>& connPICalc,
+                            const std::vector<EvalWell>& mobility,
+                            double* connPI) const
+    {
+        const auto& pu = this->phaseUsage();
+        const int   np = this->number_of_phases_;
+        for (int p = 0; p < np; ++p) {
+            // Note: E100's notion of PI value phase mobility includes
+            // the reciprocal FVF.
+            const auto connMob =
+                mobility[ flowPhaseToEbosCompIdx(p) ].value()
+                * fs.invB(p).value();
+
+            connPI[p] = connPICalc(connMob);
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+            FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
+        {
+            const auto io = pu.phase_pos[Oil];
+            const auto ig = pu.phase_pos[Gas];
+
+            const auto vapoil = connPI[ig] * fs.Rv().value();
+            const auto disgas = connPI[io] * fs.Rs().value();
+
+            connPI[io] += vapoil;
+            connPI[ig] += disgas;
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    computeConnLevelInjInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
+                           const Phase preferred_phase,
+                           const std::function<double(const double)>& connIICalc,
+                           const std::vector<EvalWell>& mobility,
+                           double* connII,
+                           DeferredLogger& deferred_logger) const
+    {
+        // Assumes single phase injection
+        const auto& pu = this->phaseUsage();
+
+        auto phase_pos = 0;
+        if (preferred_phase == Phase::GAS) {
+            phase_pos = pu.phase_pos[Gas];
+        }
+        else if (preferred_phase == Phase::OIL) {
+            phase_pos = pu.phase_pos[Oil];
+        }
+        else if (preferred_phase == Phase::WATER) {
+            phase_pos = pu.phase_pos[Water];
+        }
+        else {
+            OPM_DEFLOG_THROW(Opm::NotImplemented,
+                             "Unsupported Injector Type ("
+                             << static_cast<int>(preferred_phase)
+                             << ") for well " << this->name()
+                             << " during connection I.I. calculation",
+                             deferred_logger);
+        }
+
+        const auto zero   = EvalWell { 0.0 };
+        const auto mt     = std::accumulate(mobility.begin(), mobility.end(), zero);
+        connII[phase_pos] = connIICalc(mt.value() * fs.invB(phase_pos).value());
+    }
 } // namespace Opm

opm/simulators/wells/StandardWell.hpp

@@ -27,9 +27,10 @@
 #include <opm/simulators/linalg/bda/WellContributions.hpp>
 #endif
 
+#include <opm/simulators/wells/GasLiftRuntime.hpp>
 #include <opm/simulators/wells/RateConverter.hpp>
 #include <opm/simulators/wells/WellInterface.hpp>
-#include <opm/simulators/wells/GasLiftRuntime.hpp>
+#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 
 #include <opm/models/blackoil/blackoilpolymermodules.hh>
 #include <opm/models/blackoil/blackoilsolventmodules.hh>
@@ -38,6 +39,7 @@
 #include <opm/models/blackoil/blackoilbrinemodules.hh>
 
 #include <opm/material/densead/DynamicEvaluation.hpp>
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/ScheduleTypes.hpp>
 
 #include <dune/common/dynvector.hh>
@@ -224,6 +226,11 @@ namespace Opm
                                                  const WellState& well_state,
                                                  Opm::DeferredLogger& deferred_logger) override; // should be const?
 
+        virtual void updateProductivityIndex(const Simulator& ebosSimulator,
+                                             const WellProdIndexCalculator& wellPICalc,
+                                             WellState& well_state,
+                                             DeferredLogger& deferred_logger) const override;
+
         virtual void  addWellContributions(SparseMatrixAdapter& mat) const override;
 
         // iterate well equations with the specified control until converged
@@ -304,8 +311,19 @@ namespace Opm
         virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
                                                             DeferredLogger& deferred_logger) const override;
 
-    protected:
+        void computeConnLevelProdInd(const FluidState& fs,
+                                     const std::function<double(const double)>& connPICalc,
+                                     const std::vector<EvalWell>& mobility,
+                                     double* connPI) const;
 
+        void computeConnLevelInjInd(const typename StandardWell<TypeTag>::FluidState& fs,
+                                    const Phase preferred_phase,
+                                    const std::function<double(const double)>& connIICalc,
+                                    const std::vector<EvalWell>& mobility,
+                                    double* connII,
+                                    DeferredLogger& deferred_logger) const;
+
+    protected:
         // protected functions from the Base class
         using Base::getAllowCrossFlow;
         using Base::flowPhaseToEbosCompIdx;
@@ -315,7 +333,6 @@ namespace Opm
         using Base::wpolymer;
         using Base::wfoam;
         using Base::scalingFactor;
-        using Base::scaleProductivityIndex;
         using Base::mostStrictBhpFromBhpLimits;
 
         // protected member variables from the Base class

opm/simulators/wells/StandardWell_impl.hpp

@@ -24,6 +24,10 @@
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/linalg/MatrixBlock.hpp>
 
+#include <algorithm>
+#include <functional>
+#include <numeric>
+
 namespace Opm
 {
 
@@ -602,9 +606,6 @@ namespace Opm
         well_state.wellDissolvedGasRates()[index_of_well_] = 0.;
 
         const int np = number_of_phases_;
-        for (int p = 0; p < np; ++p) {
-            well_state.productivityIndex()[np*index_of_well_ + p] = 0.;
-        }
 
         std::vector<RateVector> connectionRates = connectionRates_; // Copy to get right size.
         for (int perf = 0; perf < number_of_perforations_; ++perf) {
@@ -683,8 +684,6 @@ namespace Opm
         } catch( ... ) {
             OPM_DEFLOG_THROW(Opm::NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
         }
-
-
     }
 
 
@@ -703,7 +702,6 @@ namespace Opm
                         Opm::DeferredLogger& deferred_logger) const
     {
         const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator);
-        const int np = number_of_phases_;
         const EvalWell& bhp = getBhp();
         const int cell_idx = well_cells_[perf];
         const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
@@ -860,27 +858,6 @@ namespace Opm
 
         // Store the perforation pressure for later usage.
         well_state.perfPress()[first_perf_ + perf] = well_state.bhp()[index_of_well_] + perf_pressure_diffs_[perf];
-
-        // Compute Productivity index if asked for
-        const auto& pu = phaseUsage();
-        const Opm::SummaryConfig& summaryConfig = ebosSimulator.vanguard().summaryConfig();
-        const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule();
-        for (int p = 0; p < np; ++p) {
-            if ( (pu.phase_pos[Water] == p && (summaryConfig.hasSummaryKey("WPIW:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name())))
-                 || (pu.phase_pos[Oil] == p && (summaryConfig.hasSummaryKey("WPIO:" + name()) || summaryConfig.hasSummaryKey("WPIL:" + name())))
-                 || (pu.phase_pos[Gas] == p && summaryConfig.hasSummaryKey("WPIG:" + name()))) {
-
-                const unsigned int compIdx = flowPhaseToEbosCompIdx(p);
-                const auto& fs = intQuants.fluidState();
-                Eval perf_pressure = getPerfCellPressure(fs);
-                const double drawdown  = well_state.perfPress()[first_perf_ + perf] - perf_pressure.value();
-                const bool new_well = schedule.hasWellGroupEvent(name(), ScheduleEvents::NEW_WELL, current_step_);
-                double productivity_index = cq_s[compIdx].value() / drawdown;
-                scaleProductivityIndex(perf, productivity_index, new_well, deferred_logger);
-                well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index;
-            }
-        }
-
     }
 
 
@@ -2307,6 +2284,80 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateProductivityIndex(const Simulator& ebosSimulator,
+                            const WellProdIndexCalculator& wellPICalc,
+                            WellState& well_state,
+                            DeferredLogger& deferred_logger) const
+    {
+        auto fluidState = [&ebosSimulator, this](const int perf)
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            return ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+        };
+
+        const int np = this->number_of_phases_;
+        auto setToZero = [np](double* x) -> void
+        {
+            std::fill_n(x, np, 0.0);
+        };
+
+        auto addVector = [np](const double* src, double* dest) -> void
+        {
+            std::transform(src, src + np, dest, dest, std::plus<>{});
+        };
+
+        auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
+        auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
+
+        setToZero(wellPI);
+
+        const auto preferred_phase = this->well_ecl_.getPreferredPhase();
+
+        const auto& allConn = this->well_ecl_.getConnections();
+        const auto  nPerf   = allConn.size();
+        auto subsetPerfID   = 0*nPerf;
+        for (auto allPerfID = 0*nPerf; allPerfID < nPerf; ++allPerfID) {
+            if (allConn[allPerfID].state() == Connection::State::SHUT) {
+                continue;
+            }
+
+            auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
+            {
+                return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
+            };
+
+            std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.0});
+            getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
+
+            const auto& fs = fluidState(subsetPerfID);
+            setToZero(connPI);
+
+            if (this->isInjector()) {
+                this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
+                                             mob, connPI, deferred_logger);
+            }
+            else {  // Production or zero flow rate
+                this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
+            }
+
+            addVector(connPI, wellPI);
+
+            ++subsetPerfID;
+            connPI += np;
+        }
+
+        assert (static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
+                "Internal logic error in processing connections for PI/II");
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
@@ -4210,4 +4261,80 @@ namespace Opm
     }
 
 
+
+
+
+    template <typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    computeConnLevelProdInd(const typename StandardWell<TypeTag>::FluidState& fs,
+                            const std::function<double(const double)>& connPICalc,
+                            const std::vector<EvalWell>& mobility,
+                            double* connPI) const
+    {
+        const auto& pu = this->phaseUsage();
+        const int   np = this->number_of_phases_;
+        for (int p = 0; p < np; ++p) {
+            // Note: E100's notion of PI value phase mobility includes
+            // the reciprocal FVF.
+            const auto connMob =
+                mobility[ flowPhaseToEbosCompIdx(p) ].value() * fs.invB(p).value();
+
+            connPI[p] = connPICalc(connMob);
+        }
+
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
+            FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
+        {
+            const auto io = pu.phase_pos[Oil];
+            const auto ig = pu.phase_pos[Gas];
+
+            const auto vapoil = connPI[ig] * fs.Rv().value();
+            const auto disgas = connPI[io] * fs.Rs().value();
+
+            connPI[io] += vapoil;
+            connPI[ig] += disgas;
+        }
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    computeConnLevelInjInd(const typename StandardWell<TypeTag>::FluidState& fs,
+                           const Phase preferred_phase,
+                           const std::function<double(const double)>& connIICalc,
+                           const std::vector<EvalWell>& mobility,
+                           double* connII,
+                           DeferredLogger& deferred_logger) const
+    {
+        // Assumes single phase injection
+        const auto& pu = this->phaseUsage();
+
+        auto phase_pos = 0;
+        if (preferred_phase == Phase::GAS) {
+            phase_pos = pu.phase_pos[Gas];
+        }
+        else if (preferred_phase == Phase::OIL) {
+            phase_pos = pu.phase_pos[Oil];
+        }
+        else if (preferred_phase == Phase::WATER) {
+            phase_pos = pu.phase_pos[Water];
+        }
+        else {
+            OPM_DEFLOG_THROW(Opm::NotImplemented,
+                             "Unsupported Injector Type ("
+                             << static_cast<int>(preferred_phase)
+                             << ") for well " << this->name()
+                             << " during connection I.I. calculation",
+                             deferred_logger);
+        }
+
+        const auto zero   = EvalWell { this->numWellEq_ + this->numEq, 0.0 };
+        const auto mt     = std::accumulate(mobility.begin(), mobility.end(), zero);
+        connII[phase_pos] = connIICalc(mt.value() * fs.invB(phase_pos).value());
+    }
 } // namespace Opm

opm/simulators/wells/WellInterface.hpp

@@ -39,6 +39,7 @@
 #include <opm/simulators/wells/VFPProperties.hpp>
 #include <opm/simulators/wells/WellHelpers.hpp>
 #include <opm/simulators/wells/WellGroupHelpers.hpp>
+#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
 #include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
 
@@ -221,6 +222,11 @@ namespace Opm
                                                  const WellState& well_state,
                                                  Opm::DeferredLogger& deferred_logger) = 0; // should be const?
 
+        virtual void updateProductivityIndex(const Simulator& ebosSimulator,
+                                             const WellProdIndexCalculator& wellPICalc,
+                                             WellState& well_state,
+                                             DeferredLogger& deferred_logger) const = 0;
+
         /// \brief Wether the Jacobian will also have well contributions in it.
         virtual bool jacobianContainsWellContributions() const
         {
@@ -513,14 +519,8 @@ namespace Opm
                                  const std::vector<double>& B_avg,
                                  Opm::DeferredLogger& deferred_logger);
 
-        void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) const;
-
         void initCompletions();
 
-        // count the number of times an output log message is created in the productivity
-        // index calculations
-        mutable int well_productivity_index_logger_counter_;
-
         bool checkConstraints(WellState& well_state,
                               const Schedule& schedule,
                               const SummaryState& summaryState,

opm/simulators/wells/WellInterface_impl.hpp

@@ -85,8 +85,6 @@ namespace Opm
 
         connectionRates_.resize(number_of_perforations_);
 
-        well_productivity_index_logger_counter_ = 0;
-
         wellIsStopped_ = false;
         if (well.getStatus() == Well::Status::STOP) {
             wellIsStopped_ = true;
@@ -1380,40 +1378,6 @@ namespace Opm
         }
     }
 
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger) const
-    {
-        const auto& connection = well_ecl_.getConnections()[(*perf_data_)[perfIdx].ecl_index];
-        if (well_ecl_.getDrainageRadius() < 0) {
-            if (new_well && perfIdx == 0) {
-                deferred_logger.warning("PRODUCTIVITY_INDEX_WARNING", "Negative drainage radius not supported. The productivity index is set to zero");
-            }
-            productivity_index = 0.0;
-            return;
-        }
-
-        if (connection.r0() > well_ecl_.getDrainageRadius()) {
-            if (new_well && well_productivity_index_logger_counter_ < 1) {
-                deferred_logger.info("PRODUCTIVITY_INDEX_INFO", "The effective radius is larger than the well drainage radius for well " + name() +
-                             " They are set to equal in the well productivity index calculations");
-                well_productivity_index_logger_counter_++;
-            }
-            return;
-        }
-
-        // For zero drainage radius the productivity index is just the transmissibility times the mobility
-        if (well_ecl_.getDrainageRadius() == 0) {
-            return;
-        }
-
-        // Scale the productivity index to account for the drainage radius.
-        // Assumes steady radial flow only valied for horizontal wells
-        productivity_index *=
-        (std::log(connection.r0() / connection.rw()) + connection.skinFactor()) /
-        (std::log(well_ecl_.getDrainageRadius() / connection.rw()) + connection.skinFactor());
-    }
-
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::addCellRates(RateVector& rates, int cellIdx) const

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -29,15 +29,15 @@
 
 #include <opm/common/ErrorMacros.hpp>
 
-#include <vector>
-#include <cassert>
-#include <string>
-#include <utility>
-#include <map>
 #include <algorithm>
 #include <array>
+#include <cassert>
 #include <iostream>
+#include <map>
 #include <numeric>
+#include <string>
+#include <utility>
+#include <vector>
 
 namespace Opm
 {
@@ -161,6 +161,7 @@ namespace Opm
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
             productivity_index_.resize(nw * np, 0.0);
+            conn_productivity_index_.resize(nperf * np, 0.0);
             well_potentials_.resize(nw * np, 0.0);
 
             perfRatePolymer_.clear();
@@ -515,16 +516,20 @@ namespace Opm
 
             using rt = data::Rates::opt;
             std::vector< rt > phs( np );
+            std::vector<rt> pi(np);
             if( pu.phase_used[Water] ) {
                 phs.at( pu.phase_pos[Water] ) = rt::wat;
+                pi .at( pu.phase_pos[Water] ) = rt::productivity_index_water;
             }
 
             if( pu.phase_used[Oil] ) {
                 phs.at( pu.phase_pos[Oil] ) = rt::oil;
+                pi .at( pu.phase_pos[Oil] ) = rt::productivity_index_oil;
             }
 
             if( pu.phase_used[Gas] ) {
                 phs.at( pu.phase_pos[Gas] ) = rt::gas;
+                pi .at( pu.phase_pos[Gas] ) = rt::productivity_index_gas;
             }
 
             /* this is a reference or example on **how** to convert from
@@ -612,13 +617,14 @@ namespace Opm
 
                 size_t local_comp_index = 0;
                 for( auto& comp : well.connections) {
-                    const auto rates = this->perfPhaseRates().begin()
-                                     + (np * wt.second[ 1 ])
-                                     + (np * local_comp_index);
+                    const auto connPhaseOffset = np * (wt.second[1] + local_comp_index);
 
+                    const auto rates  = this->perfPhaseRates().begin() + connPhaseOffset;
+                    const auto connPI = this->connectionProductivityIndex().begin() + connPhaseOffset;
 
                     for( int i = 0; i < np; ++i ) {
-                        comp.rates.set( phs[ i ], *(rates + i) );
+                        comp.rates.set( phs[ i ], *(rates  + i) );
+                        comp.rates.set( pi [ i ], *(connPI + i) );
                     }
                     if ( pu.has_polymer ) {
                         comp.rates.set( rt::polymer, this->perfRatePolymer()[wt.second[1] + local_comp_index]);
@@ -981,6 +987,14 @@ namespace Opm
             return productivity_index_;
         }
 
+        std::vector<double>& connectionProductivityIndex() {
+            return this->conn_productivity_index_;
+        }
+
+        const std::vector<double>& connectionProductivityIndex() const {
+            return this->conn_productivity_index_;
+        }
+
         std::vector<double>& wellPotentials() {
             return well_potentials_;
         }
@@ -1262,6 +1276,9 @@ namespace Opm
         // Productivity Index
         std::vector<double> productivity_index_;
 
+        // Connection-level Productivity Index
+        std::vector<double> conn_productivity_index_;
+
         // Well potentials
         std::vector<double> well_potentials_;