Merge pull request #6114 from atgeirr/safer_well_state

Store primary variables in well state, and improve well initialization

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -429,16 +429,14 @@ namespace Opm {
         // we need the inj_multiplier from the previous time step
         this->initInjMult();
 
-        const auto& summaryState = simulator_.vanguard().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.
             // This is done only for producers, as injectors will only have a single
             // nonzero phase anyway.
             for (const auto& well : well_container_) {
-                const bool zero_target = well->stoppedOrZeroRateTarget(simulator_, this->wellState(), local_deferredLogger);
-                if (well->isProducer() && !zero_target) {
-                    well->updateWellStateRates(simulator_, this->wellState(), local_deferredLogger);
+                if (well->isProducer() && !well->wellIsStopped()) {
+                    well->initializeProducerWellState(simulator_, this->wellState(), local_deferredLogger);
                 }
             }
         }
@@ -466,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,
@@ -559,7 +557,7 @@ namespace Opm {
 
             // initialize rates/previous rates to prevent zero fractions in vfp-interpolation
             if (well->isProducer()) {
-                well->updateWellStateRates(simulator_, 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

@@ -111,7 +111,16 @@ update(const WellState<Scalar>& well_state,
         if (stop_or_zero_rate_target && seg == 0) {
             value_[seg][WQTotal] = 0;
         }
-        if (std::abs(total_seg_rate) > 0.) {
+
+        // todo to map old fraction to new perforations for now start from scratch.
+        if (ws.primaryvar.size() == value_.size() * numWellEq) {
+            if (has_wfrac_variable) {
+                value_[seg][WFrac] = ws.primaryvar[seg * numWellEq + WFrac];
+            }
+            if (has_gfrac_variable) {
+                value_[seg][GFrac] = ws.primaryvar[seg * numWellEq + GFrac];
+            }
+        } else if (std::abs(total_seg_rate) > 0.) {
             if (has_wfrac_variable) {
                 const int water_pos = pu.phase_pos[Water];
                 value_[seg][WFrac] = well_.scalingFactor(water_pos) * segment_rates[well_.numPhases() * seg + water_pos] / total_seg_rate;
@@ -230,6 +239,15 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
     const int oil_pos = pu.phase_pos[Oil];
 
     auto& ws = well_state.well(well_.indexOfWell());
+
+    // Store primary variables
+    ws.primaryvar.resize(value_.size() * numWellEq);
+    for (std::size_t seg = 0; seg < value_.size(); ++seg) {
+        for (int ii = 0; ii < numWellEq; ++ii) {
+            ws.primaryvar[seg * numWellEq + ii] = value_[seg][ii];
+        }
+    }
+
     auto& segments = ws.segments;
     auto& segment_rates = segments.rates;
     auto& disgas = segments.dissolved_gas_rate;

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1483,6 +1483,8 @@ namespace Opm
                 return false;
         }
 
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+
         std::vector<std::vector<Scalar> > residual_history;
         std::vector<Scalar> measure_history;
         int it = 0;
@@ -1601,6 +1603,8 @@ namespace Opm
                 return false;
         }
 
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+
         std::vector<std::vector<Scalar> > residual_history;
         std::vector<Scalar> measure_history;
         int it = 0;

opm/simulators/wells/SingleWellState.cpp

@@ -252,6 +252,7 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
         return;
     }
 
+    std::fill(this->surface_rates.begin(), this->surface_rates.end(), 0.0);
     switch (prod_controls.cmode) {
     case Well::ProducerCMode::ORAT:
         assert(this->pu.phase_used[BlackoilPhases::Liquid]);
@@ -268,15 +269,12 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
     case Well::ProducerCMode::GRUP:
     case Well::ProducerCMode::THP:
     case Well::ProducerCMode::BHP:
-        if (this->pu.phase_used[BlackoilPhases::Liquid]) {
-            this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
-        }
-        if (this->pu.phase_used[BlackoilPhases::Aqua]) {
-            this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
-        }
-        if (this->pu.phase_used[BlackoilPhases::Vapour]){
-            this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -100000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
-        }
+        // Keeping all rates at zero, they will be initialized properly in
+        // 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::initializeProducerWellState().
         break;
 
     default:

opm/simulators/wells/SingleWellState.hpp

@@ -65,7 +65,7 @@ public:
         serializer(well_potentials);
         serializer(productivity_index);
         serializer(implicit_ipr_a);
-        serializer(implicit_ipr_b);        
+        serializer(implicit_ipr_b);
         serializer(surface_rates);
         serializer(reservoir_rates);
         serializer(prev_surface_rates);
@@ -76,6 +76,7 @@ public:
         serializer(production_cmode);
         serializer(filtrate_conc);
         serializer(perf_data);
+        serializer(primaryvar);
     }
 
     bool operator==(const SingleWellState&) const;
@@ -115,7 +116,7 @@ public:
     Events events;
     WellInjectorCMode injection_cmode{WellInjectorCMode::CMODE_UNDEFINED};
     WellProducerCMode production_cmode{WellProducerCMode::CMODE_UNDEFINED};
-
+    std::vector<Scalar> primaryvar;
 
     /// Special purpose method to support dynamically rescaling a well's
     /// CTFs through WELPI.

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -163,7 +163,17 @@ update(const WellState<Scalar>& well_state,
             }
     }
 
-    if (std::abs(total_well_rate) > 0.) {
+    if (ws.primaryvar.size() > 0) {
+        if constexpr (has_wfrac_variable) {
+             value_[WFrac] = ws.primaryvar[WFrac];
+        }
+        if constexpr (has_gfrac_variable) {
+             value_[GFrac] = ws.primaryvar[GFrac];
+        }
+        if constexpr (Indices::enableSolvent) {
+            value_[SFrac] = ws.primaryvar[SFrac];
+         }
+    } else if (std::abs(total_well_rate) > 0.) {
         if constexpr (has_wfrac_variable) {
             value_[WFrac] = well_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
         }
@@ -175,7 +185,6 @@ update(const WellState<Scalar>& well_state,
         if constexpr (Indices::enableSolvent) {
             value_[SFrac] = well_.scalingFactor(Indices::contiSolventEqIdx) * ws.sum_solvent_rates() / total_well_rate ;
         }
-
     } else { // total_well_rate == 0
         if (well_.isInjector()) {
             // only single phase injection handled
@@ -399,6 +408,7 @@ copyToWellState(WellState<Scalar>& well_state,
     }
 
     auto& ws = well_state.well(well_.indexOfWell());
+    ws.primaryvar = value_;
     ws.bhp = value_[Bhp];
 
     // calculate the phase rates based on the primary variables

opm/simulators/wells/StandardWell_impl.hpp

@@ -2324,6 +2324,8 @@ namespace Opm
                              const GroupState<Scalar>& group_state,
                              DeferredLogger& deferred_logger)
     {
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+
         const int max_iter = this->param_.max_inner_iter_wells_;
         int it = 0;
         bool converged;
@@ -2371,6 +2373,8 @@ namespace Opm
                                const bool fixed_control /*false*/,
                                const bool fixed_status /*false*/)
     {
+        updatePrimaryVariables(simulator, well_state, deferred_logger);
+
         const int max_iter = this->param_.max_inner_iter_wells_;
         int it = 0;
         bool converged = false;

opm/simulators/wells/WellInterface.hpp

@@ -327,10 +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 updateWellStateRates(const Simulator& simulator,
-                              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

@@ -1137,6 +1137,10 @@ namespace Opm
         const auto& summaryState = simulator.vanguard().summaryState();
         const auto& schedule = simulator.vanguard().schedule();
 
+        // Discard old primary variables, the new well state
+        // may not be anywhere near the old one.
+        ws.primaryvar.resize(0);
+
         if (this->wellIsStopped()) {
             for (int p = 0; p<np; ++p) {
                 ws.surface_rates[p] = 0;
@@ -1603,10 +1607,11 @@ namespace Opm
     template <typename TypeTag>
     void
     WellInterface<TypeTag>::
-    updateWellStateRates(const Simulator& simulator,
-                         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.
@@ -1624,29 +1629,65 @@ namespace Opm
             }
         }
 
-        // Calculate the rates that follow from the current primary variables.
-        std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
+        // 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
+            // No nonzero rates on input.
+            // 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 Scalar 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.
+        // If we are here, we had a single nonzero rate for the well,
+        // typically from a rate constraint. We must make sure it is
+        // respected, so if it was lower than the calculated rate for
+        // the same phase we scale all rates to match.
         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) {
+        const Scalar computed_rate = well_q_s[nonzero_rate_index];
+        if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) {
+            // Note that both rates below are negative. The factor should be < 1.0.
+            const Scalar factor = initial_nonzero_rate / computed_rate;
+            assert(factor < 1.0);
             for (int p = 0; p < this->number_of_phases_; ++p) {
+                // We skip the nonzero_rate_index, as that should remain as it was.
                 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]);
+                    ws.surface_rates[p] = factor * well_q_s[p];
                 }
             }
+            return;
+        }
+
+        // If we are here, we had a single nonzero rate, but it was
+        // higher than the one calculated from the bhp limit, so we
+        // use the calculated rates.
+        for (int p = 0; p < this->number_of_phases_; ++p) {
+            ws.surface_rates[p] = 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;