Unconditionally Calculate PI at End of Timestep

This commit ensures that we calculate the well and connection level
per-phase steady-state productivity index (PI) at the end of a
completed time step (triggered from endTimeStep()).

We add a new data member,

    BlackoilWellModel<>::prod_index_calc_

which holds one WellProdIndexCalculator for each of the process'
local wells and a new interface member function

    WellInterface::updateProductivityIndex

which uses a per-well PI calculator to actually compute the PI
values and store those in the WellState.  Implement this member
function for both StandardWell and MultisegmentWell.  Were it not
for 'getMobility' existing only in the derived classes, the two
equal implementations could be merged and moved to the interface.

We also add a new data member to the WellStateFullyImplicitBlackoil
to hold the connection-level PI values.  Finally, remove the
conditional PI calculation from StandardWell's well equation
assembly routine.

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

@@ -44,7 +44,6 @@ namespace Opm
         using typename Base::RateConverterType;
         using typename Base::SparseMatrixAdapter;
 
-
         /// the number of reservior equations
         using Base::numEq;
 
@@ -168,6 +167,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

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1176,6 +1176,94 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    updateProductivityIndex(const Simulator& ebosSimulator,
+                            const WellProdIndexCalculator& wellPICalc,
+                            WellState& well_state,
+                            DeferredLogger&) const
+    {
+        auto recipFVF = [&ebosSimulator, this](const int perf, const int p) -> double
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            const auto& fs = ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+
+            return fs.invB(p).value();
+        };
+
+        auto Rs = [&ebosSimulator, this](const int perf) -> double
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            const auto& fs = ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+
+            return fs.Rs().value();
+        };
+
+        auto Rv = [&ebosSimulator, this](const int perf) -> double
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            const auto& fs = ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+
+            return fs.Rv().value();
+        };
+
+        const auto& pu = this->phaseUsage();
+        const int   np = this->number_of_phases_;
+
+        auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
+        auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
+
+        for (int p = 0; p < np; ++p) {
+            wellPI[p] = 0.0;
+        }
+
+        const auto& allConn = this->well_ecl_.getConnections();
+        const auto  nPerf   = allConn.size();
+        for (auto allPerfID = 0*nPerf, subsetPerfID = 0*nPerf; allPerfID < nPerf; ++allPerfID) {
+            if (allConn[allPerfID].state() == Connection::State::SHUT) {
+                continue;
+            }
+
+            std::vector<EvalWell> mob(num_components_, 0.0);
+            getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob);
+
+            for (int p = 0; p < np; ++p) {
+                // Note: E100's notion of PI value phase mobility includes
+                // the reciprocal FVF.
+                const auto connMob = mob[ flowPhaseToEbosCompIdx(p) ].value() * recipFVF(subsetPerfID, p);
+                connPI[p] = wellPICalc.connectionProdIndStandard(allPerfID, 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] * Rv(subsetPerfID);
+                const auto disgas = connPI[io] * Rs(subsetPerfID);
+
+                connPI[io] += vapoil;
+                connPI[ig] += disgas;
+            }
+
+            for (int p = 0; p < np; ++p) {
+                wellPI[p] += connPI[p];
+            }
+
+            ++subsetPerfID;
+            connPI += np;
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::

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>
@@ -224,6 +225,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
@@ -315,7 +321,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

@@ -602,9 +602,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 +680,6 @@ namespace Opm
         } catch( ... ) {
             OPM_DEFLOG_THROW(Opm::NumericalIssue,"Error when inverting local well equations for well " + name(), deferred_logger);
         }
-
-
     }
 
 
@@ -703,7 +698,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 +854,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 +2280,94 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateProductivityIndex(const Simulator& ebosSimulator,
+                            const WellProdIndexCalculator& wellPICalc,
+                            WellState& well_state,
+                            DeferredLogger& deferred_logger) const
+    {
+        auto recipFVF = [&ebosSimulator, this](const int perf, const int p) -> double
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            const auto& fs = ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+
+            return fs.invB(p).value();
+        };
+
+        auto Rs = [&ebosSimulator, this](const int perf) -> double
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            const auto& fs = ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+
+            return fs.Rs().value();
+        };
+
+        auto Rv = [&ebosSimulator, this](const int perf) -> double
+        {
+            const auto cell_idx = this->well_cells_[perf];
+            const auto& fs = ebosSimulator.model()
+               .cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
+
+            return fs.Rv().value();
+        };
+
+        const auto& pu = this->phaseUsage();
+        const int   np = this->number_of_phases_;
+
+        auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
+        auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
+
+        for (int p = 0; p < np; ++p) {
+            wellPI[p] = 0.0;
+        }
+
+        const auto& allConn = this->well_ecl_.getConnections();
+        const auto  nPerf   = allConn.size();
+        for (auto allPerfID = 0*nPerf, subsetPerfID = 0*nPerf; allPerfID < nPerf; ++allPerfID) {
+            if (allConn[allPerfID].state() == Connection::State::SHUT) {
+                continue;
+            }
+
+            std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.0});
+            getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
+
+            for (int p = 0; p < np; ++p) {
+                // Note: E100's notion of PI value phase mobility includes
+                // the reciprocal FVF.
+                const auto connMob = mob[ flowPhaseToEbosCompIdx(p) ].value() * recipFVF(subsetPerfID, p);
+                connPI[p] = wellPICalc.connectionProdIndStandard(allPerfID, 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] * Rv(subsetPerfID);
+                const auto disgas = connPI[io] * Rs(subsetPerfID);
+
+                connPI[io] += vapoil;
+                connPI[ig] += disgas;
+            }
+
+            for (int p = 0; p < np; ++p) {
+                wellPI[p] += connPI[p];
+            }
+
+            ++subsetPerfID;
+            connPI += np;
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::

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_;