Refine producer initialization.

Instead of evaluating the well at the state given by its current
state (primary variables) use the bhp limit (or 1 bar) to get the
rates that provide the correct fractions.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -429,7 +429,6 @@ namespace Opm {
         // we need the inj_multiplier from the previous time step
         this->initInjMult();
 
-        const auto& summaryState = this->summaryState();
         if (alternative_well_rate_init_) {
             // Update the well rates of well_state_, if only single-phase rates, to
             // have proper multi-phase rates proportional to rates at bhp zero.
@@ -437,7 +436,7 @@ namespace Opm {
             // nonzero phase anyway.
             for (const auto& well : well_container_) {
                 if (well->isProducer()) {
-                    well->initializeWellState(simulator_, this->groupState(), summaryState, this->wellState(), local_deferredLogger);
+                    well->initializeProducerWellState(simulator_, this->wellState(), local_deferredLogger);
                 }
             }
         }
@@ -465,7 +464,7 @@ namespace Opm {
         const Group& fieldGroup = this->schedule().getGroup("FIELD", reportStepIdx);
         WellGroupHelpers<Scalar>::updateGuideRates(fieldGroup,
                                                    this->schedule(),
-                                                   summaryState,
+                                                   this->summaryState(),
                                                    this->phase_usage_,
                                                    reportStepIdx,
                                                    simulationTime,
@@ -558,7 +557,7 @@ namespace Opm {
 
             // initialize rates/previous rates to prevent zero fractions in vfp-interpolation
             if (well->isProducer()) {
-                well->initializeWellState(simulator_, this->groupState(), this->summaryState(), this->wellState(), deferred_logger);
+                well->initializeProducerWellState(simulator_, this->wellState(), deferred_logger);
             }
             if (well->isVFPActive(deferred_logger)) {
                 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -112,9 +112,6 @@ update(const WellState<Scalar>& well_state,
             value_[seg][WQTotal] = 0;
         }
 
-        // Production wells should have been initialized using reservoir data.
-        assert(well_.isProducer() == ws.initialized_from_reservoir);
-
         // todo to map old fraction to new perforations for now start from scratch.
         if (ws.primaryvar.size() == value_.size() * numWellEq) {
             if (has_wfrac_variable) {

opm/simulators/wells/SingleWellState.cpp

@@ -51,7 +51,6 @@ SingleWellState(const std::string& name_,
     , prev_surface_rates(pu_.num_phases)
     , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases)
     , trivial_group_target(false)
-    , initialized_from_reservoir(false)
 {
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
@@ -264,11 +263,11 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
     case Well::ProducerCMode::THP:
     case Well::ProducerCMode::BHP:
         // Setting all rates to zero, they will be initialized properly in
-        // a call to WellInterface::initializeWellState() later, which will
+        // a call to WellInterface::initializeProducerWellState() later, which will
         // use the reservoir state to find a better initial value.
         // This also applies to the ORAT/WRAT/GRAT above, but then only the
         // rate not set in the above will be modified in
-        // WellInterface::initializeWellState().
+        // WellInterface::initializeProducerWellState().
         if (this->pu.phase_used[BlackoilPhases::Liquid]) {
             this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = 0.0;
         }

opm/simulators/wells/SingleWellState.hpp

@@ -77,7 +77,6 @@ public:
         serializer(filtrate_conc);
         serializer(perf_data);
         serializer(primaryvar);
-        serializer(initialized_from_reservoir);
     }
 
     bool operator==(const SingleWellState&) const;
@@ -118,7 +117,6 @@ public:
     WellInjectorCMode injection_cmode{WellInjectorCMode::CMODE_UNDEFINED};
     WellProducerCMode production_cmode{WellProducerCMode::CMODE_UNDEFINED};
     std::vector<Scalar> primaryvar;
-    bool initialized_from_reservoir{false};
 
     /// Special purpose method to support dynamically rescaling a well's
     /// CTFs through WELPI.

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -163,9 +163,6 @@ update(const WellState<Scalar>& well_state,
             }
     }
 
-    // Production wells should have been initialized using reservoir data.
-    assert(well_.isProducer() == ws.initialized_from_reservoir);
-
     if (ws.primaryvar.size() > 0) {
         if constexpr (has_wfrac_variable) {
              value_[WFrac] = ws.primaryvar[WFrac];

opm/simulators/wells/WellInterface.hpp

@@ -327,12 +327,10 @@ public:
 
     /// Modify the well_state's rates if there is only one nonzero rate.
     /// If so, that rate is kept as is, but the others are set proportionally
-    /// to the rates returned by computeCurrentWellRates().
-    void initializeWellState(const Simulator& simulator,
-                             const GroupState<Scalar>& group_state,
-                             const SummaryState& summary_state,
-                             WellState<Scalar>& well_state,
-                             DeferredLogger& deferred_logger) const;
+    /// to the rates at bhp limit or bhp 1 bar.
+    void initializeProducerWellState(const Simulator& simulator,
+                                     WellState<Scalar>& well_state,
+                                     DeferredLogger& deferred_logger) const;
 
     void solveWellEquation(const Simulator& simulator,
                            WellState<Scalar>& well_state,

opm/simulators/wells/WellInterface_impl.hpp

@@ -1603,12 +1603,11 @@ namespace Opm
     template <typename TypeTag>
     void
     WellInterface<TypeTag>::
-    initializeWellState(const Simulator& simulator,
-                        const GroupState<Scalar>& /*group_state*/,
-                        const SummaryState& /* summary_state */,
-                        WellState<Scalar>& well_state,
-                        DeferredLogger& deferred_logger) const
+    initializeProducerWellState(const Simulator& simulator,
+                                WellState<Scalar>& well_state,
+                                DeferredLogger& deferred_logger) const
     {
+        assert(this->isProducer());
         OPM_TIMEFUNCTION();
         // Check if the rates of this well only are single-phase, do nothing
         // if more than one nonzero rate.
@@ -1626,29 +1625,48 @@ namespace Opm
             }
         }
 
-        // Calculate the rates that follow from the current primary variables.
-        // TODO switch to ccomputeWellRatesWithBhp().
-        std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
-        ws.initialized_from_reservoir = true;
+        // Calculate rates at bhp limit, or 1 bar if no limit.
+        std::vector<Scalar> well_q_s(this->number_of_phases_, 0.0);
+        bool rates_evaluated_at_1bar = false;
+        {
+            const auto& summary_state = simulator.vanguard().summaryState();
+            const auto& prod_controls = this->well_ecl_.productionControls(summary_state);
+            const double bhp_limit = std::max(prod_controls.bhp_limit, 1.0 * unit::barsa);
+            this->computeWellRatesWithBhp(simulator, bhp_limit, well_q_s, deferred_logger);
+            // Remember of we evaluated the rates at (approx.) 1 bar or not.
+            rates_evaluated_at_1bar = (bhp_limit < 1.1 * unit::barsa);
+            // Check that no rates are positive.
+            if (std::any_of(well_q_s.begin(), well_q_s.end(), [](Scalar q) { return q > 0.0; })) {
+                // Did we evaluate at 1 bar? If not, then we can try again at 1 bar.
+                if (!rates_evaluated_at_1bar) {
+                    this->computeWellRatesWithBhp(simulator, 1.0 * unit::barsa, well_q_s, deferred_logger);
+                    rates_evaluated_at_1bar = true;
+                }
+                // At this point we can only set the wrong-direction (if any) values to zero.
+                for (auto& q : well_q_s) {
+                    q = std::min(q, Scalar{0.0});
+                }
+            }
+        }
 
         if (nonzero_rate_index == -1) {
             // No nonzero rates.
-            // Use the computed rate directly
+            // Use the computed rate directly, or scaled by a factor
+            // 0.5 (to avoid too high values) if it was evaluated at 1 bar.
+            const double factor = rates_evaluated_at_1bar ? 0.5 : 1.0;
             for (int p = 0; p < this->number_of_phases_; ++p) {
-               ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
+               ws.surface_rates[p] = factor * well_q_s[p];
             }
             return;
         }
 
         // Set the currently-zero phase flows to be nonzero in proportion to well_q_s.
         const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index];
-        const int comp_idx_nz = this->flowPhaseToModelCompIdx(nonzero_rate_index);
-        if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) {
+        if (std::abs(well_q_s[nonzero_rate_index]) > floating_point_error_epsilon) {
             for (int p = 0; p < this->number_of_phases_; ++p) {
                 if (p != nonzero_rate_index) {
-                    const int comp_idx = this->flowPhaseToModelCompIdx(p);
                     Scalar& rate = ws.surface_rates[p];
-                    rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
+                    rate = (initial_nonzero_rate / well_q_s[nonzero_rate_index]) * (well_q_s[p]);
                 }
             }
         }

opm/simulators/wells/WellState.cpp

@@ -409,6 +409,7 @@ void WellState<Scalar>::init(const std::vector<Scalar>& cellPressures,
             }
 
             new_well.surface_rates = prev_well.surface_rates;
+            new_well.prev_surface_rates = prev_well.prev_surface_rates;
             new_well.reservoir_rates = prev_well.reservoir_rates;
             new_well.well_potentials = prev_well.well_potentials;