Revise Mixture Density Method for No-Flow Producers

This commit switches the approach introduced in commit eeb1b7e36 (PR
#3169) to using a mobility weighted average of cell level densities
for the connection level mixture densities in no-flow producing
wells.  We also use the recent stoppedOrZeroRateTarget() predicate
to identify those no-flow producing wells instead of inspecting the
connection flow rates.

The mobility weighted average gives a more monotone pressure buildup
for the stopped wells and this is usually what the engineer wants.
This revised approach furthermore needs fewer cell-level dynamic
properties so simplify the computeProperties() signature by
introducing a structure for the property callback functions and
update the callers accordingly.

opm/simulators/wells/StandardWellConnections.cpp

@@ -419,6 +419,44 @@ initialiseConnectionMixture(const int                  num_comp,
     }
 }
 
+template <class FluidSystem, class Indices>
+void StandardWellConnections<FluidSystem, Indices>::
+computeDensitiesForStoppedProducer(const DensityPropertyFunctions& prop_func)
+{
+    const auto np = this->well_.numPhases();
+
+    const auto modPhIx = [this, np]() {
+        auto phIx = std::vector<int>(np);
+
+        for (auto p = 0*np; p < np; ++p) {
+            phIx[p] = this->well_.flowPhaseToModelPhaseIdx(p);
+        }
+
+        return phIx;
+    }();
+
+    auto mob = std::vector<Scalar>(np);
+    auto rho = std::vector<Scalar>(np);
+
+    const auto nperf = this->well_.numPerfs();
+    for (auto perf = 0*nperf; perf < nperf; ++perf) {
+        const auto cell_idx = this->well_.cells()[perf];
+
+        prop_func.mobility     (cell_idx, modPhIx, mob);
+        prop_func.densityInCell(cell_idx, modPhIx, rho);
+
+        auto& rho_i = this->perf_densities_[perf];
+
+        auto mt = rho_i = Scalar{0};
+        for (auto p = 0*np; p < np; ++p) {
+            rho_i += mob[p] * rho[p];
+            mt    += mob[p];
+        }
+
+        rho_i /= mt;
+    }
+}
+
 template<class FluidSystem, class Indices>
 typename StandardWellConnections<FluidSystem, Indices>::Properties
 StandardWellConnections<FluidSystem,Indices>::
@@ -625,85 +663,27 @@ copyInPerforationRates(const Properties&       props,
 
 template<class FluidSystem, class Indices>
 void StandardWellConnections<FluidSystem,Indices>::
-computeProperties(const WellState<Scalar>& well_state,
-                  const std::function<Scalar(int,int)>& invB,
-                  const std::function<Scalar(int,int)>& mobility,
-                  const std::function<Scalar(int)>& solventInverseFormationVolumeFactor,
-                  const std::function<Scalar(int)>& solventMobility,
-                  const Properties& props,
-                  DeferredLogger& deferred_logger)
+computeProperties(const bool                      stopped_or_zero_rate_target,
+                  const WellState<Scalar>&        well_state,
+                  const DensityPropertyFunctions& prop_func,
+                  const Properties&               props,
+                  DeferredLogger&                 deferred_logger)
 {
-    // Compute densities
-    const int nperf = well_.numPerfs();
-    const int np = well_.numPhases();
-    const auto& ws = well_state.well(well_.indexOfWell());
-    const auto& perf_data = ws.perf_data;
-    const auto& perf_rates_state = perf_data.phase_rates;
+    // Step 1: Compute mixture densities in well-bore.  Result values stored
+    // in data member perf_densities_.
+    if (stopped_or_zero_rate_target && this->well_.isProducer()) {
+        this->computeDensitiesForStoppedProducer(prop_func);
+    }
+    else {
+        // Injector or flowing producer.
+        const auto perfRates = this->copyInPerforationRates
+            (props, well_state.well(this->well_.indexOfWell()).perf_data);
 
-    auto perfRates = this->copyInPerforationRates
-        (props, well_state.well(this->well_.indexOfWell()).perf_data);
-
-    // for producers where all perforations have zero rate we
-    // approximate the perforation mixture using the mobility ratio
-    // and weight the perforations using the well transmissibility.
-    bool all_zero = std::all_of(perfRates.begin(), perfRates.end(),
-                                [](Scalar val) { return val == 0.0; });
-    const auto& comm = well_.parallelWellInfo().communication();
-    if (comm.size() > 1)
-    {
-        all_zero =  (comm.min(all_zero ? 1 : 0) == 1);
+        this->computeDensities(perfRates, props, deferred_logger);
     }
 
-    if (all_zero && well_.isProducer()) {
-        Scalar total_tw = 0;
-        for (int perf = 0; perf < nperf; ++perf) {
-            total_tw += well_.wellIndex()[perf];
-        }
-        if (comm.size() > 1)
-        {
-            total_tw = comm.sum(total_tw);
-        }
-        // for producers where all perforations have zero rates we
-        // approximate the perforation mixture ration using the (invB * mobility) ratio,
-        // and weight the perforation rates using the well transmissibility.
-        for (int perf = 0; perf < nperf; ++perf) {
-            const int cell_idx = well_.cells()[perf];
-            const Scalar well_tw_fraction = well_.wellIndex()[perf] / total_tw;
-            Scalar total_mobility = 0.0;
-            Scalar total_invB = 0.;
-            for (int p = 0; p < np; ++p) {
-                int modelPhaseIdx = well_.flowPhaseToModelPhaseIdx(p);
-                total_mobility += invB(cell_idx, modelPhaseIdx) * mobility(cell_idx, modelPhaseIdx);
-                total_invB += invB(cell_idx, modelPhaseIdx);
-            }
-            if constexpr (Indices::enableSolvent) {
-                total_mobility += solventInverseFormationVolumeFactor(cell_idx) * solventMobility(cell_idx);
-                total_invB += solventInverseFormationVolumeFactor(cell_idx);
-            }
-            const bool non_zero_total_mobility = total_mobility > std::numeric_limits<Scalar>::min();
-            assert(total_invB > std::numeric_limits<Scalar>::min());
-            // for the perforation having zero mobility for all the phases, we use a small value to generate a small
-            // perforation rates for those perforations, at the same time, we can use the rates to recover the mixing
-            // ratios for those perforations.
-            constexpr Scalar small_value = 1.e-10;
-            for (int p = 0; p < np; ++p) {
-                const int modelPhaseIdx = well_.flowPhaseToModelPhaseIdx(p);
-                const auto mob_ratio = non_zero_total_mobility
-                                       ? mobility(cell_idx, modelPhaseIdx) / total_mobility
-                                       : small_value / total_invB;
-                perfRates[perf * well_.numComponents() + p] = well_tw_fraction * invB(cell_idx, modelPhaseIdx) * mob_ratio;
-            }
-            if constexpr (Indices::enableSolvent) {
-                const auto mob_ratio = non_zero_total_mobility
-                                        ? solventMobility(cell_idx) / total_mobility
-                                        : small_value / total_invB;
-                perfRates[perf * well_.numComponents() + Indices::contiSolventEqIdx] =
-                        well_tw_fraction * solventInverseFormationVolumeFactor(cell_idx) * mob_ratio;
-            }
-        }
-    }
-
-    this->computeDensities(perfRates, props, deferred_logger);
+    // Step 2: Use those mixture densities to calculate pressure drops along
+    // well-bore.  Result values stored in data member perf_pressure_diffs_.
     this->computePressureDelta();
 }
 

opm/simulators/wells/StandardWellConnections.hpp

@@ -66,18 +66,22 @@ public:
         std::function<Scalar(int)>     solventRefDensity{};
     };
 
+    struct DensityPropertyFunctions
+    {
+        std::function<void(int, const std::vector<int>&, std::vector<Scalar>&)> mobility{};
+        std::function<void(int, const std::vector<int>&, std::vector<Scalar>&)> densityInCell{};
+    };
+
     Properties
     computePropertiesForPressures(const WellState<Scalar>&         well_state,
                                   const PressurePropertyFunctions& propFunc) const;
 
     //! \brief Compute connection properties (densities, pressure drop, ...)
-    void computeProperties(const WellState<Scalar>& well_state,
-                           const std::function<Scalar(int,int)>& invB,
-                           const std::function<Scalar(int,int)>& mobility,
-                           const std::function<Scalar(int)>& solventInverseFormationVolumeFactor,
-                           const std::function<Scalar(int)>& solventMobility,
-                           const Properties& props,
-                           DeferredLogger& deferred_logger);
+    void computeProperties(const bool                      stop_or_zero_rate_target,
+                           const WellState<Scalar>&        well_state,
+                           const DensityPropertyFunctions& prop_func,
+                           const Properties&               props,
+                           DeferredLogger&                 deferred_logger);
 
     //! \brief Returns density for first perforation.
     Scalar rho() const
@@ -140,6 +144,8 @@ private:
                           const Properties& props,
                           DeferredLogger& deferred_logger);
 
+    void computeDensitiesForStoppedProducer(const DensityPropertyFunctions& prop_func);
+
     std::vector<Scalar>
     calculatePerforationOutflow(const std::vector<Scalar>& perfComponentRates) const;
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -1304,41 +1304,49 @@ namespace Opm
 
 
     template<typename TypeTag>
-    void
-    StandardWell<TypeTag>::
+    void StandardWell<TypeTag>::
     computeWellConnectionDensitesPressures(const Simulator& simulator,
                                            const WellState<Scalar>& well_state,
                                            const WellConnectionProps& props,
                                            DeferredLogger& deferred_logger)
     {
-        std::function<Scalar(int,int)> invB =
-        [&simulator](int cell_idx, int phase_idx)
-        {
-            return simulator.model().intensiveQuantities(cell_idx, 0).fluidState().invB(phase_idx).value();
-        };
-        std::function<Scalar(int,int)> mobility =
-        [&simulator](int cell_idx, int phase_idx)
-        {
-            return simulator.model().intensiveQuantities(cell_idx, 0).mobility(phase_idx).value();
-        };
-        std::function<Scalar(int)> invFac =
-        [&simulator](int cell_idx)
-        {
-            return simulator.model().intensiveQuantities(cell_idx, 0).solventInverseFormationVolumeFactor().value();
-        };
-        std::function<Scalar(int)> solventMobility =
-        [&simulator](int cell_idx)
-        {
-            return simulator.model().intensiveQuantities(cell_idx, 0).solventMobility().value();
+        // Cell level dynamic property call-back functions as fall-back
+        // option for calculating connection level mixture densities in
+        // stopped or zero-rate producer wells.
+        const auto prop_func = typename StdWellEval::StdWellConnections::DensityPropertyFunctions {
+            // This becomes slightly more palatable with C++20's designated
+            // initialisers.
+
+            // mobility: Phase mobilities in specified cell.
+            [&model = simulator.model()](const int               cell,
+                                         const std::vector<int>& phases,
+                                         std::vector<Scalar>&    mob)
+            {
+                const auto& iq = model.intensiveQuantities(cell, /* time_idx = */ 0);
+
+                std::transform(phases.begin(), phases.end(), mob.begin(),
+                               [&iq](const int phase) { return iq.mobility(phase).value(); });
+            },
+
+            // densityInCell: Reservoir condition phase densities in
+            // specified cell.
+            [&model = simulator.model()](const int               cell,
+                                         const std::vector<int>& phases,
+                                         std::vector<Scalar>&    rho)
+            {
+                const auto& fs = model.intensiveQuantities(cell, /* time_idx = */ 0).fluidState();
+
+                std::transform(phases.begin(), phases.end(), rho.begin(),
+                               [&fs](const int phase) { return fs.density(phase).value(); });
+            }
         };
 
-        this->connections_.computeProperties(well_state,
-                                             invB,
-                                             mobility,
-                                             invFac,
-                                             solventMobility,
-                                             props,
-                                             deferred_logger);
+        const auto stopped_or_zero_rate_target = this->
+            stoppedOrZeroRateTarget(simulator, well_state, deferred_logger);
+
+        this->connections_
+            .computeProperties(stopped_or_zero_rate_target, well_state,
+                               prop_func, props, deferred_logger);
     }