Merge pull request #4694 from bska/wbp-compute-values

Hook New WBPn Calculation Up to Well Model

ebos/eclwriter.hh

@@ -210,7 +210,7 @@ public:
             simulator_.vanguard().setupTime();
 
         const auto localWellData            = simulator_.problem().wellModel().wellData();
-        const auto localWBP                 = data::WellBlockAveragePressures{};
+        const auto localWBP                 = simulator_.problem().wellModel().wellBlockAveragePressures();
         const auto localGroupAndNetworkData = simulator_.problem().wellModel()
             .groupAndNetworkData(reportStepNum);
 

opm/simulators/wells/BlackoilWellModel.hpp

@@ -38,35 +38,39 @@
 #include <unordered_map>
 #include <vector>
 
+#include <opm/input/eclipse/Schedule/Group/Group.hpp>
+#include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
 #include <opm/input/eclipse/Schedule/Schedule.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
-#include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
-#include <opm/input/eclipse/Schedule/Group/Group.hpp>
 
-#include <opm/simulators/timestepping/SimulatorReport.hpp>
 #include <opm/simulators/flow/countGlobalCells.hpp>
 #include <opm/simulators/flow/SubDomain.hpp>
+
 #include <opm/simulators/wells/BlackoilWellModelGeneric.hpp>
 #include <opm/simulators/wells/BlackoilWellModelGuideRates.hpp>
+#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
 #include <opm/simulators/wells/GasLiftSingleWell.hpp>
-#include <opm/simulators/wells/GasLiftWellState.hpp>
 #include <opm/simulators/wells/GasLiftSingleWellGeneric.hpp>
 #include <opm/simulators/wells/GasLiftStage2.hpp>
-#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
+#include <opm/simulators/wells/GasLiftWellState.hpp>
+#include <opm/simulators/wells/MultisegmentWell.hpp>
+#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
+#include <opm/simulators/wells/ParallelWellInfo.hpp>
 #include <opm/simulators/wells/PerforationData.hpp>
-#include <opm/simulators/wells/VFPInjProperties.hpp>
-#include <opm/simulators/wells/VFPProdProperties.hpp>
-#include <opm/simulators/wells/WellState.hpp>
-#include <opm/simulators/wells/WGState.hpp>
 #include <opm/simulators/wells/RateConverter.hpp>
 #include <opm/simulators/wells/RegionAverageCalculator.hpp>
-#include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/wells/StandardWell.hpp>
-#include <opm/simulators/wells/MultisegmentWell.hpp>
+#include <opm/simulators/wells/VFPInjProperties.hpp>
+#include <opm/simulators/wells/VFPProdProperties.hpp>
+#include <opm/simulators/wells/WGState.hpp>
 #include <opm/simulators/wells/WellGroupHelpers.hpp>
+#include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
-#include <opm/simulators/wells/ParallelWellInfo.hpp>
+#include <opm/simulators/wells/WellState.hpp>
+
+#include <opm/simulators/timestepping/SimulatorReport.hpp>
 #include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
+
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmatrix.hh>
@@ -266,6 +270,11 @@ namespace Opm {
                 return wsrpt;
             }
 
+            data::WellBlockAveragePressures wellBlockAveragePressures() const
+            {
+                return this->computeWellBlockAveragePressures();
+            }
+
             // subtract Binv(D)rw from r;
             void apply( BVector& r) const;
 
@@ -387,6 +396,13 @@ namespace Opm {
             std::unique_ptr<RateConverterType> rateConverter_{};
             std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};
 
+            struct WBPCalcID
+            {
+                std::optional<typename std::vector<WellInterfacePtr>::size_type> openWellIdx_{};
+                std::size_t wbpCalcIdx_{};
+            };
+
+            std::vector<WBPCalcID> wbpCalcMap_{};
 
             SimulatorReportSingle last_report_{};
 
@@ -447,6 +463,16 @@ namespace Opm {
             // setting the well_solutions_ based on well_state.
             void updatePrimaryVariables(DeferredLogger& deferred_logger);
 
+            void initializeWBPCalculationService();
+
+            data::WellBlockAveragePressures
+            computeWellBlockAveragePressures() const;
+
+            ParallelWBPCalculation::EvaluatorFactory
+            makeWellSourceEvaluatorFactory(const std::vector<Well>::size_type wellIdx) const;
+
+            void registerOpenWellsForWBPCalculation();
+
             void updateAverageFormationFactor();
 
             void computePotentials(const std::size_t widx,

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -49,6 +49,7 @@
 #include <opm/simulators/wells/BlackoilWellModelGuideRates.hpp>
 #include <opm/simulators/wells/BlackoilWellModelRestart.hpp>
 #include <opm/simulators/wells/GasLiftStage2.hpp>
+#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
 #include <opm/simulators/wells/VFPProperties.hpp>
 #include <opm/simulators/wells/WellFilterCake.hpp>
 #include <opm/simulators/wells/WellGroupHelpers.hpp>
@@ -81,6 +82,7 @@ BlackoilWellModelGeneric(Schedule& schedule,
     , eclState_(eclState)
     , comm_(comm)
     , phase_usage_(phase_usage)
+    , wbpCalculationService_ { eclState.gridDims(), comm_ }
     , guideRate_(schedule)
     , active_wgstate_(phase_usage)
     , last_valid_wgstate_(phase_usage)
@@ -287,57 +289,89 @@ BlackoilWellModelGeneric::
 initializeWellPerfData()
 {
     well_perf_data_.resize(wells_ecl_.size());
+
+    this->conn_idx_map_.clear();
+    this->conn_idx_map_.reserve(wells_ecl_.size());
+
     int well_index = 0;
     for (const auto& well : wells_ecl_) {
         int connection_index = 0;
+
         // INVALID_ECL_INDEX marks no above perf available
         int connection_index_above = ParallelWellInfo::INVALID_ECL_INDEX;
+
         well_perf_data_[well_index].clear();
         well_perf_data_[well_index].reserve(well.getConnections().size());
-        CheckDistributedWellConnections checker(well, local_parallel_well_info_[well_index].get());
+
+        auto& connIdxMap = this->conn_idx_map_
+            .emplace_back(well.getConnections().size());
+
+        CheckDistributedWellConnections checker {
+            well, this->local_parallel_well_info_[well_index].get()
+        };
+
         bool hasFirstConnection = false;
         bool firstOpenConnection = true;
+
         auto& parallelWellInfo = this->local_parallel_well_info_[well_index].get();
         parallelWellInfo.beginReset();
 
         for (const auto& connection : well.getConnections()) {
-            const int active_index = compressedIndexForInterior(connection.global_index());
+            const auto active_index =
-            if (connection.state() == Connection::State::OPEN) {
+                this->compressedIndexForInterior(connection.global_index());
+
+            const auto connIsOpen =
+                connection.state() == Connection::State::OPEN;
+
+            if (active_index >= 0) {
+                connIdxMap.addActiveConnection(connection_index, connIsOpen);
+            }
+
+            if ((connIsOpen && (active_index >= 0)) || !connIsOpen) {
+                checker.connectionFound(connection_index);
+            }
+
+            if (connIsOpen) {
                 if (active_index >= 0) {
-                    if (firstOpenConnection)
+                    if (firstOpenConnection) {
-                    {
                         hasFirstConnection = true;
                     }
-                    checker.connectionFound(connection_index);
+
-                    PerforationData pd;
+                    auto pd = PerforationData{};
                     pd.cell_index = active_index;
                     pd.connection_transmissibility_factor = connection.CF();
                     pd.satnum_id = connection.satTableId();
                     pd.ecl_index = connection_index;
+
                     well_perf_data_[well_index].push_back(pd);
+
                     parallelWellInfo.pushBackEclIndex(connection_index_above,
                                                       connection_index);
                 }
+
                 firstOpenConnection = false;
-                // Next time this index is the one above as each open connection is
+
-                // is stored somehwere.
+                // Next time this index is the one above as each open
+                // connection is stored somewhere.
                 connection_index_above = connection_index;
-            } else {
-                checker.connectionFound(connection_index);
-                if (connection.state() != Connection::State::SHUT) {
-                    OPM_THROW(std::runtime_error,
-                              "Connection state: " +
-                              Connection::State2String(connection.state()) +
-                              " not handled");
-                }
             }
-            // Note: we rely on the connections being filtered! I.e. there are only connections
+            else if (connection.state() != Connection::State::SHUT) {
-            // to active cells in the global grid.
+                OPM_THROW(std::runtime_error,
+                          fmt::format("Connection state '{}' not handled",
+                                      Connection::State2String(connection.state())));
+            }
+
+            // Note: we rely on the connections being filtered!  I.e., there
+            // are only connections to active cells in the global grid.
             ++connection_index;
         }
+
         parallelWellInfo.endReset();
+
         checker.checkAllConnectionsFound();
+
         parallelWellInfo.communicateFirstPerforation(hasFirstConnection);
+
         ++well_index;
     }
 }

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -25,20 +25,28 @@
 
 #include <opm/output/data/GuideRateValue.hpp>
 
-#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
 #include <opm/input/eclipse/Schedule/Group/GuideRate.hpp>
+#include <opm/input/eclipse/Schedule/Well/PAvg.hpp>
+#include <opm/input/eclipse/Schedule/Well/PAvgCalculator.hpp>
+#include <opm/input/eclipse/Schedule/Well/PAvgCalculatorCollection.hpp>
+#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
+#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
+#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
 #include <opm/simulators/wells/PerforationData.hpp>
 #include <opm/simulators/wells/WellFilterCake.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/WGState.hpp>
 
+#include <cstddef>
 #include <functional>
 #include <map>
 #include <memory>
+#include <optional>
 #include <string>
+#include <type_traits>
 #include <unordered_map>
 #include <unordered_set>
 #include <vector>
@@ -419,6 +427,94 @@ protected:
 
     std::vector<Well> wells_ecl_;
     std::vector<std::vector<PerforationData>> well_perf_data_;
+
+    /// Connection index mappings
+    class ConnectionIndexMap
+    {
+    public:
+        /// Constructor.
+        ///
+        /// \param[in] numConns Total number of well connections, both open
+        ///   and closed/shut.  Typically \code WellConnections::size() \endcode.
+        explicit ConnectionIndexMap(const std::size_t numConns)
+            : local_(numConns, -1)
+        {
+            this->global_.reserve(numConns);
+            this->open_.reserve(numConns);
+        }
+
+        /// Enumerate/map new active connection.
+        ///
+        /// \param[in] connIdx Global well connection index.  Must be an
+        ///   integer in the range 0..numConns-1.
+        ///
+        /// \param[in] connIsOpen Whether or not the connection is
+        ///   open/flowing.
+        void addActiveConnection(const int  connIdx,
+                                 const bool connIsOpen)
+        {
+            this->local_[connIdx] =
+                static_cast<int>(this->global_.size());
+
+            this->global_.push_back(connIdx);
+
+            const auto open_conn_idx = connIsOpen
+                ? this->num_open_conns_++
+                : -1;
+
+            this->open_.push_back(open_conn_idx);
+        }
+
+        /// Get local connection IDs/indices of every existing well
+        /// connection.
+        ///
+        /// Negative value (-1) for connections that don't intersect the
+        /// current rank.
+        const std::vector<int>& local() const
+        {
+            return this->local_;
+        }
+
+        /// Get global connection ID of local (on-rank) connection.
+        ///
+        /// \param[in] connIdx Local connection index.
+        ///
+        /// \return Global connection ID of \p connIdx.
+        int global(const int connIdx) const
+        {
+            return this->global_[connIdx];
+        }
+
+        /// Get open connection ID of local (on-rank) connection.
+        ///
+        /// \param[in] connIdx Local connection index.
+        ///
+        /// \return Open connection ID of \p connIdx.  Integer in the range
+        ///   0..#open connections - 1 if the connection is open or negative
+        ///   value (-1) otherwise.
+        int open(const int connIdx) const
+        {
+            return this->open_[connIdx];
+        }
+
+    private:
+        /// Local connection IDs/indices of every existing well connection.
+        /// Negative value (-1) for connections that don't intersect the
+        /// current rank.
+        std::vector<int> local_{};
+
+        /// Global connection index of each on-rank reservoir connection.
+        /// Reverse/transpose mapping of \c local_.
+        std::vector<int> global_{};
+
+        /// Open connection index of each on-rank reservoir connection.
+        std::vector<int> open_{};
+
+        /// Number of open connections on this rank.
+        int num_open_conns_{0};
+    };
+
+    std::vector<ConnectionIndexMap> conn_idx_map_{};
     std::function<bool(const Well&)> not_on_process_{};
 
     // a vector of all the wells.
@@ -430,6 +526,7 @@ protected:
     std::vector<std::reference_wrapper<ParallelWellInfo>> local_parallel_well_info_;
 
     std::vector<WellProdIndexCalculator> prod_index_calc_;
+    mutable ParallelWBPCalculation wbpCalculationService_;
 
     std::vector<int> pvt_region_idx_;
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -27,8 +27,11 @@
 #include <opm/input/eclipse/Units/UnitSystem.hpp>
 #include <opm/input/eclipse/Schedule/Network/Balance.hpp>
 #include <opm/input/eclipse/Schedule/Network/ExtNetwork.hpp>
+#include <opm/input/eclipse/Schedule/Well/PAvgDynamicSourceData.hpp>
 
 #include <opm/simulators/wells/BlackoilWellModelConstraints.hpp>
+#include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
+#include <opm/simulators/wells/ParallelWBPCalculation.hpp>
 #include <opm/simulators/wells/VFPProperties.hpp>
 #include <opm/simulators/utils/MPIPacker.hpp>
 #include <opm/simulators/linalg/bda/WellContributions.hpp>
@@ -73,6 +76,48 @@ namespace Opm {
 
         this->alternative_well_rate_init_ =
             EWOMS_GET_PARAM(TypeTag, bool, AlternativeWellRateInit);
+
+        this->wbpCalculationService_
+            .localCellIndex([this](const std::size_t globalIndex)
+            { return this->compressedIndexForInterior(globalIndex); })
+            .evalCellSource([this](const int                                     localCell,
+                                   PAvgDynamicSourceData::SourceDataSpan<double> sourceTerms)
+            {
+                using Item = PAvgDynamicSourceData::SourceDataSpan<double>::Item;
+
+                const auto* intQuants = this->ebosSimulator_.model()
+                    .cachedIntensiveQuantities(localCell, /*timeIndex = */0);
+                const auto& fs = intQuants->fluidState();
+
+                sourceTerms.set(Item::PoreVol, intQuants->porosity().value() *
+                                this->ebosSimulator_.model().dofTotalVolume(localCell));
+
+                constexpr auto io = FluidSystem::oilPhaseIdx;
+                constexpr auto ig = FluidSystem::gasPhaseIdx;
+                constexpr auto iw = FluidSystem::waterPhaseIdx;
+
+                // Ideally, these would be 'constexpr'.
+                const auto haveOil = FluidSystem::phaseIsActive(io);
+                const auto haveGas = FluidSystem::phaseIsActive(ig);
+                const auto haveWat = FluidSystem::phaseIsActive(iw);
+
+                auto weightedPhaseDensity = [&fs](const auto ip)
+                {
+                    return fs.saturation(ip).value() * fs.density(ip).value();
+                };
+
+                if (haveOil)      { sourceTerms.set(Item::Pressure, fs.pressure(io).value()); }
+                else if (haveGas) { sourceTerms.set(Item::Pressure, fs.pressure(ig).value()); }
+                else              { sourceTerms.set(Item::Pressure, fs.pressure(iw).value()); }
+
+                // Strictly speaking, assumes SUM(s[p]) == 1.
+                auto rho = 0.0;
+                if (haveOil) { rho += weightedPhaseDensity(io); }
+                if (haveGas) { rho += weightedPhaseDensity(ig); }
+                if (haveWat) { rho += weightedPhaseDensity(iw); }
+
+                sourceTerms.set(Item::MixtureDensity, rho);
+            });
     }
 
     template<typename TypeTag>
@@ -209,6 +254,7 @@ namespace Opm {
             // We must therefore provide it with updated cell pressures
             this->initializeWellPerfData();
             this->initializeWellState(timeStepIdx, summaryState);
+            this->initializeWBPCalculationService();
 
             // handling MS well related
             if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model
@@ -230,11 +276,16 @@ namespace Opm {
             }
 
             {
-                const auto& sched_state = this->schedule()[timeStepIdx];
                 // update VFP properties
+                const auto& sched_state = this->schedule()[timeStepIdx];
                 vfp_properties_ = std::make_unique<VFPProperties>(sched_state.vfpinj(),
                                                                   sched_state.vfpprod(),
                                                                   this->wellState());
+            }
+
+            {
+                const auto& sched_state = this->schedule()[timeStepIdx];
+
                 this->initializeWellProdIndCalculators();
                 if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) {
                     this->runWellPIScaling(timeStepIdx, local_deferredLogger);
@@ -799,6 +850,8 @@ namespace Opm {
                 }
             }
         }
+
+        this->registerOpenWellsForWBPCalculation();
     }
 
 
@@ -1801,6 +1854,154 @@ namespace Opm {
 
 
 
+
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    initializeWBPCalculationService()
+    {
+        this->wbpCalcMap_.clear();
+        this->wbpCalcMap_.resize(this->wells_ecl_.size());
+
+        this->registerOpenWellsForWBPCalculation();
+
+        auto wellID = std::size_t{0};
+        for (const auto& well : this->wells_ecl_) {
+            this->wbpCalcMap_[wellID].wbpCalcIdx_ = this->wbpCalculationService_
+                .createCalculator(well,
+                                  this->local_parallel_well_info_[wellID],
+                                  this->conn_idx_map_[wellID].local(),
+                                  this->makeWellSourceEvaluatorFactory(wellID));
+
+            ++wellID;
+        }
+
+        this->wbpCalculationService_.defineCommunication();
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    data::WellBlockAveragePressures
+    BlackoilWellModel<TypeTag>::
+    computeWellBlockAveragePressures() const
+    {
+        auto wbpResult = data::WellBlockAveragePressures{};
+
+        using Calculated = PAvgCalculator::Result::WBPMode;
+        using Output = data::WellBlockAvgPress::Quantity;
+
+        this->wbpCalculationService_.collectDynamicValues();
+
+        const auto numWells = this->wells_ecl_.size();
+        for (auto wellID = 0*numWells; wellID < numWells; ++wellID) {
+            const auto calcIdx = this->wbpCalcMap_[wellID].wbpCalcIdx_;
+            const auto& well = this->wells_ecl_[wellID];
+
+            if (! well.hasRefDepth()) {
+                // Can't perform depth correction without at least a
+                // fall-back datum depth.
+                continue;
+            }
+
+            this->wbpCalculationService_
+                .inferBlockAveragePressures(calcIdx, well.pavg(),
+                                            this->gravity_,
+                                            well.getWPaveRefDepth());
+
+            const auto& result = this->wbpCalculationService_
+                .averagePressures(calcIdx);
+
+            auto& reported = wbpResult.values[well.name()];
+
+            reported[Output::WBP]  = result.value(Calculated::WBP);
+            reported[Output::WBP4] = result.value(Calculated::WBP4);
+            reported[Output::WBP5] = result.value(Calculated::WBP5);
+            reported[Output::WBP9] = result.value(Calculated::WBP9);
+        }
+
+        return wbpResult;
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    ParallelWBPCalculation::EvaluatorFactory
+    BlackoilWellModel<TypeTag>::
+    makeWellSourceEvaluatorFactory(const std::vector<Well>::size_type wellIdx) const
+    {
+        using Span = PAvgDynamicSourceData::SourceDataSpan<double>;
+        using Item = typename Span::Item;
+
+        return [wellIdx, this]() -> ParallelWBPCalculation::Evaluator
+        {
+            if (! this->wbpCalcMap_[wellIdx].openWellIdx_.has_value()) {
+                // Well is stopped/shut.  Return evaluator for stopped wells.
+                return []([[maybe_unused]] const int connIdx, Span sourceTerm)
+                {
+                    // Well/connection is stopped/shut.  Set all items to
+                    // zero.
+
+                    sourceTerm
+                        .set(Item::Pressure      , 0.0)
+                        .set(Item::PoreVol       , 0.0)
+                        .set(Item::MixtureDensity, 0.0);
+                };
+            }
+
+            // Well is open.  Return an evaluator for open wells/open connections.
+            return [this, wellPtr = this->well_container_[*this->wbpCalcMap_[wellIdx].openWellIdx_].get()]
+                (const int connIdx, Span sourceTerm)
+            {
+                // Note: The only item which actually matters for the WBP
+                // calculation at the well reservoir connection level is the
+                // mixture density.  Set other items to zero.
+
+                const auto& connIdxMap =
+                    this->conn_idx_map_[wellPtr->indexOfWell()];
+
+                const auto rho = wellPtr->
+                    connectionDensity(connIdxMap.global(connIdx),
+                                      connIdxMap.open(connIdx));
+
+                sourceTerm
+                    .set(Item::Pressure      , 0.0)
+                    .set(Item::PoreVol       , 0.0)
+                    .set(Item::MixtureDensity, rho);
+            };
+        };
+    }
+
+
+
+
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    registerOpenWellsForWBPCalculation()
+    {
+        assert (this->wbpCalcMap_.size() == this->wells_ecl_.size());
+
+        for (auto& wbpCalc : this->wbpCalcMap_) {
+            wbpCalc.openWellIdx_.reset();
+        }
+
+        auto openWellIdx = typename std::vector<WellInterfacePtr>::size_type{0};
+        for (const auto* openWell : this->well_container_generic_) {
+            this->wbpCalcMap_[openWell->indexOfWell()].openWellIdx_ = openWellIdx++;
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/MultisegmentWell.hpp

@@ -137,6 +137,9 @@ namespace Opm
                                              WellState& well_state,
                                              DeferredLogger& deferred_logger) const override;
 
+        double connectionDensity(const int globalConnIdx,
+                                 const int openConnIdx) const override;
+
         void addWellContributions(SparseMatrixAdapter& jacobian) const override;
 
         void addWellPressureEquations(PressureMatrix& mat,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -21,8 +21,12 @@
 #include <opm/common/Exceptions.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
 
+#include <opm/input/eclipse/Schedule/MSW/Segment.hpp>
 #include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
+#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
+#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
+
 #include <opm/input/eclipse/Units/Units.hpp>
 
 #include <opm/material/densead/EvaluationFormat.hpp>
@@ -703,6 +707,28 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    double
+    MultisegmentWell<TypeTag>::
+    connectionDensity(const int globalConnIdx,
+                      [[maybe_unused]] const int openConnIdx) const
+    {
+        // Simple approximation: Mixture density at reservoir connection is
+        // mixture density at connection's segment.
+
+        const auto segNum = this->wellEcl()
+            .getConnections()[globalConnIdx].segment();
+
+        const auto segIdx = this->wellEcl()
+            .getSegments().segmentNumberToIndex(segNum);
+
+        return this->segments_.density(segIdx).value();
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::

opm/simulators/wells/StandardWell.hpp

@@ -185,6 +185,9 @@ namespace Opm
                                              WellState& well_state,
                                              DeferredLogger& deferred_logger) const override;
 
+        virtual double connectionDensity(const int globalConnIdx,
+                                         const int openConnIdx) const override;
+
         virtual void addWellContributions(SparseMatrixAdapter& mat) const override;
 
         virtual void addWellPressureEquations(PressureMatrix& mat,

opm/simulators/wells/StandardWellConnections.hpp

@@ -73,7 +73,19 @@ public:
     //! \brief Returns density for first perforation.
     Scalar rho() const
     {
-        return this->perf_densities_.empty() ? 0.0 : perf_densities_[0];
+        return this->rho(0);
+    }
+
+    //! \brief Returns density for specific perforation/connection.
+    //!
+    //! \param[in] i Connection index
+    //!
+    //! \return Mixture density at connection \p i.
+    Scalar rho(const typename std::vector<Scalar>::size_type i) const
+    {
+        return (i < this->perf_densities_.size())
+            ? this->perf_densities_[i]
+            : 0.0;
     }
 
     //! \brief Returns pressure drop for a given perforation.

opm/simulators/wells/StandardWell_impl.hpp

@@ -1570,6 +1570,23 @@ namespace Opm
 
 
 
+
+
+    template<typename TypeTag>
+    double
+    StandardWell<TypeTag>::
+    connectionDensity([[maybe_unused]] const int globalConnIdx,
+                      const int openConnIdx) const
+    {
+        return (openConnIdx < 0)
+            ? 0.0
+            : this->connections_.rho(openConnIdx);
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::

opm/simulators/wells/WellInterface.hpp

@@ -254,6 +254,9 @@ public:
                                          WellState& well_state,
                                          DeferredLogger& deferred_logger) const = 0;
 
+    virtual double connectionDensity(const int globalConnIdx,
+                                     const int openConnIdx) const = 0;
+
     /// \brief Wether the Jacobian will also have well contributions in it.
     virtual bool jacobianContainsWellContributions() const
     {