Compute well productivity index and pass it to the output

opm/autodiff/StandardWell.hpp

@@ -187,6 +187,7 @@ namespace Opm
         using Base::wellHasTHPConstraints;
         using Base::mostStrictBhpFromBhpLimits;
         using Base::scalingFactor;
+        using Base::scaleProductivityIndex;
 
         // protected member variables from the Base class
         using Base::current_step_;

opm/autodiff/StandardWell_impl.hpp

@@ -461,6 +461,10 @@ namespace Opm
         well_state.wellVaporizedOilRates()[index_of_well_] = 0.;
         well_state.wellDissolvedGasRates()[index_of_well_] = 0.;
 
+        for (int p = 0; p < np; ++p) {
+            well_state.productivityIndex()[np*index_of_well_ + p] = 0.;
+        }
+
         for (int perf = 0; perf < number_of_perforations_; ++perf) {
 
             const int cell_idx = well_cells_[perf];
@@ -601,8 +605,19 @@ 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
+            for (int p = 0; p < np; ++p) {
+                const unsigned int compIdx = flowPhaseToEbosCompIdx(p);
+                const double drawdown = well_state.perfPress()[first_perf_ + perf] - intQuants.fluidState().pressure(FluidSystem::oilPhaseIdx).value();
+                double productivity_index = cq_s[compIdx].value() / drawdown;
+                scaleProductivityIndex(perf, productivity_index);
+                well_state.productivityIndex()[np*index_of_well_ + p] += productivity_index;
             }
 
+        }
+
+
         // add vol * dF/dt + Q to the well equations;
         for (int componentIdx = 0; componentIdx < numWellConservationEq; ++componentIdx) {
             EvalWell resWell_loc = (wellSurfaceVolumeFraction(componentIdx) - F0_[componentIdx]) * volume / dt;

opm/autodiff/WellInterface.hpp

@@ -333,6 +333,8 @@ namespace Opm
 
         void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector<double>& B_avg, bool terminal_output);
 
+        void scaleProductivityIndex(const int perfIdx, double& productivity_index) const;
+
     };
 
 }

opm/autodiff/WellInterface_impl.hpp

@@ -103,6 +103,7 @@ namespace Opm
                       saturation_table_number_.begin() );
         }
         well_efficiency_factor_ = 1.0;
+
     }
 
 
@@ -1048,4 +1049,37 @@ namespace Opm
             }
         }
     }
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::scaleProductivityIndex(const int perfIdx, double& productivity_index) const
+    {
+
+        const auto& connection = well_ecl_->getConnections(current_step_)[perfIdx];
+
+        if (well_ecl_->getDrainageRadius(current_step_) < 0) {
+            OpmLog::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(current_step_)) {
+            OpmLog::info("PRODUCTIVITY_INDEX_INFO", "The effective radius is larger then the well drainage radius for well " + name() +
+                         " They are set to equal in the well productivity index calculations");
+            return;
+        }
+
+        // For zero drainage radius the productivity index is just the transmissibility times the mobility
+        if (well_ecl_->getDrainageRadius(current_step_) == 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(current_step_) / connection.rw()) + connection.skinFactor());
+    }
+
+
 }

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -128,6 +128,8 @@ namespace Opm
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
 
+            productivity_index_.resize(nw * np, 0.0);
+
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if(prevState && !prevState->wellMap().empty()) {
@@ -262,6 +264,8 @@ namespace Opm
             well_dissolved_gas_rates_.resize(nw, 0.0);
             well_vaporized_oil_rates_.resize(nw, 0.0);
 
+            productivity_index_.resize(nw * np, 0.0);
+
             // Ensure that we start out with zero rates by default.
             perfphaserates_.clear();
             perfphaserates_.resize(nperf * np, 0.0);
@@ -480,6 +484,18 @@ namespace Opm
                     well.rates.set( rt::reservoir_gas, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Gas]] );
                 }
 
+                if ( pu.phase_used[Water] ) {
+                    well.rates.set( rt::productivity_index_water, this->productivity_index_[well_rate_index + pu.phase_pos[Water]] );
+                }
+
+                if ( pu.phase_used[Oil] ) {
+                    well.rates.set( rt::productivity_index_oil, this->productivity_index_[well_rate_index + pu.phase_pos[Oil]] );
+                }
+
+                if ( pu.phase_used[Gas] ) {
+                    well.rates.set( rt::productivity_index_gas, this->productivity_index_[well_rate_index + pu.phase_pos[Gas]] );
+                }
+
                 well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] );
                 well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] );
 
@@ -761,6 +777,14 @@ namespace Opm
             return top_segment_index_[w];
         }
 
+        std::vector<double>& productivityIndex() {
+            return productivity_index_;
+        }
+
+        const std::vector<double>& productivityIndex() const {
+            return productivity_index_;
+        }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
@@ -793,6 +817,9 @@ namespace Opm
         std::vector<int> top_segment_index_;
         int nseg_; // total number of the segments
 
+        // Productivity Index
+        std::vector<double> productivity_index_;
+
     };
 
 } // namespace Opm