Updates to clean up.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -436,8 +436,7 @@ namespace Opm {
             // 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) {
+                if (well->isProducer()) {
                     well->initializeWellState(simulator_, this->groupState(), summaryState, this->wellState(), local_deferredLogger);
                 }
             }

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -111,29 +111,29 @@ update(const WellState<Scalar>& well_state,
         if (stop_or_zero_rate_target && seg == 0) {
             value_[seg][WQTotal] = 0;
         }
-        assert(ws.initializedFromReservoir());
-        // tot to map old fraction to new perforations for now start from scratch.
-        bool prim_set = ws.multiseg_primaryvar.size() == value_.size();
-        if(prim_set){
-        //if (std::abs(total_seg_rate) > 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) {
-                value_[seg][WFrac] = ws.multiseg_primaryvar[seg][WFrac];
+                value_[seg][WFrac] = ws.primaryvar[seg * numWellEq + WFrac];
             }
             if (has_gfrac_variable) {
-                value_[seg][GFrac] = ws.multiseg_primaryvar[seg][GFrac];
+                value_[seg][GFrac] = ws.primaryvar[seg * numWellEq + GFrac];
             }
-        } else {
-        if (std::abs(total_seg_rate) > 0.) {
+        } 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;
             }
             if (has_gfrac_variable) {
                 const int gas_pos = pu.phase_pos[Gas];
-                value_[seg][GFrac] = well_.scalingFactor(gas_pos) * segment_rates[well_.numPhases() * seg + gas_pos] / total_seg_rate; 
+                value_[seg][GFrac] = well_.scalingFactor(gas_pos) * segment_rates[well_.numPhases() * seg + gas_pos] / total_seg_rate;
             }
             // what about water and gas injection?
-         } else { // total_seg_rate == 0
+        } else { // total_seg_rate == 0
             if (well_.isInjector()) {
                 // only single phase injection handled
                 auto phase = well.getInjectionProperties().injectorType;
@@ -164,7 +164,6 @@ update(const WellState<Scalar>& well_state,
                 }
             }
         }
-        }
     }
     setEvaluationsFromValues();
 }
@@ -237,16 +236,22 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
     static constexpr int Gas = BlackoilPhases::Vapour;
     static constexpr int Oil = BlackoilPhases::Liquid;
     static constexpr int Water = BlackoilPhases::Aqua;
-    
+
     const auto pvtReg = std::max(well_.wellEcl().pvt_table_number() - 1, 0);
 
     const PhaseUsage& pu = well_.phaseUsage();
     const int oil_pos = pu.phase_pos[Oil];
 
     auto& ws = well_state.well(well_.indexOfWell());
-    if constexpr ( numWellEq == 4) {
-        ws.multiseg_primaryvar = value_;
+
+    // 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/SingleWellState.cpp

@@ -51,7 +51,7 @@ 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)
+    , 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;
@@ -263,14 +263,20 @@ update_producer_targets(const Well& ecl_well, const SummaryState& st)
     case Well::ProducerCMode::GRUP:
     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
+        // 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().
         if (this->pu.phase_used[BlackoilPhases::Liquid]) {
-            this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = 0;//-1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = 0.0;
         }
         if (this->pu.phase_used[BlackoilPhases::Aqua]) {
-            this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = 0;//-1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = 0.0;
         }
         if (this->pu.phase_used[BlackoilPhases::Vapour]){
-            this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = 0;//-100000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = 0.0;
         }
         break;
 

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,9 +76,8 @@ public:
         serializer(production_cmode);
         serializer(filtrate_conc);
         serializer(perf_data);
-        serializer(stw_primaryvar);
-        serializer(multiseg_primaryvar);
-        serializer(initialized_from_reservoir_);
+        serializer(primaryvar);
+        serializer(initialized_from_reservoir);
     }
 
     bool operator==(const SingleWellState&) const;
@@ -118,7 +117,8 @@ public:
     Events events;
     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.
@@ -153,13 +153,8 @@ public:
 
     Scalar sum_filtrate_rate() const;
     Scalar sum_filtrate_total() const;
-    std::vector<Scalar> stw_primaryvar;
-    std::vector<std::array<Scalar, 4>> multiseg_primaryvar;
-    bool initializedFromReservoir() const { return initialized_from_reservoir_; }
-   void setInitializedFromReservoir(bool value) { initialized_from_reservoir_ = value; }
-private:
-    bool initialized_from_reservoir_ = false;
 
+private:
     Scalar sum_connection_rates(const std::vector<Scalar>& connection_rates) const;
 };
 

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -162,30 +162,28 @@ update(const WellState<Scalar>& well_state,
                 value_[WQTotal] = 0.;
             }
     }
-    assert(ws.initializedFromReservoir());
-    //if (std::abs(total_well_rate) > 0.) {
-        
-    bool prim_set = ws.stw_primaryvar.size()>0;
-    if(prim_set){   
+
+    // 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.stw_primaryvar[WFrac];
+             value_[WFrac] = ws.primaryvar[WFrac];
         }
         if constexpr (has_gfrac_variable) {
-             value_[GFrac] = ws.stw_primaryvar[GFrac];
+             value_[GFrac] = ws.primaryvar[GFrac];
         }
         if constexpr (Indices::enableSolvent) {
-            value_[SFrac] = ws.stw_primaryvar[SFrac];
+            value_[SFrac] = ws.primaryvar[SFrac];
          }
-     }else{
-     assert(prim_set==false);   
-     if (std::abs(total_well_rate) > 0.) {
+    } 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;
         }
         if constexpr (has_gfrac_variable) {
             value_[GFrac] = well_.scalingFactor(pu.phase_pos[Gas]) *
-                             (ws.surface_rates[pu.phase_pos[Gas]] -
-                              (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
+                            (ws.surface_rates[pu.phase_pos[Gas]] -
+                             (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0) ) / total_well_rate ;
         }
         if constexpr (Indices::enableSolvent) {
             value_[SFrac] = well_.scalingFactor(Indices::contiSolventEqIdx) * ws.sum_solvent_rates() / total_well_rate ;
@@ -234,8 +232,8 @@ update(const WellState<Scalar>& well_state,
             OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well", deferred_logger);
         }
     }
-    }
-    //BHP
+
+    // BHP
     value_[Bhp] = ws.bhp;
     setEvaluationsFromValues();
 }
@@ -348,6 +346,7 @@ copyToWellState(WellState<Scalar>& well_state,
     static constexpr int Water = BlackoilPhases::Aqua;
     static constexpr int Oil = BlackoilPhases::Liquid;
     static constexpr int Gas = BlackoilPhases::Vapour;
+
     const PhaseUsage& pu = well_.phaseUsage();
     std::vector<Scalar> F(well_.numPhases(), 0.0);
     [[maybe_unused]] Scalar F_solvent = 0.0;
@@ -410,9 +409,9 @@ copyToWellState(WellState<Scalar>& well_state,
         F_solvent /= well_.scalingFactor(Indices::contiSolventEqIdx);
         F[pu.phase_pos[Gas]] += F_solvent;
     }
-    
+
     auto& ws = well_state.well(well_.indexOfWell());
-    ws.stw_primaryvar = value_;
+    ws.primaryvar = value_;
     ws.bhp = value_[Bhp];
 
     // calculate the phase rates based on the primary variables

opm/simulators/wells/WellInterface.hpp

@@ -334,10 +334,6 @@ public:
                              WellState<Scalar>& well_state,
                              DeferredLogger& deferred_logger) const;
 
-    // void updateWellStateRates(const Simulator& simulator,
-    //                           WellState<Scalar>& well_state,
-    //                           DeferredLogger& deferred_logger) const;
-
     void solveWellEquation(const Simulator& simulator,
                            WellState<Scalar>& well_state,
                            const GroupState<Scalar>& group_state,

opm/simulators/wells/WellInterface_impl.hpp

@@ -1598,57 +1598,6 @@ namespace Opm
         return scaling_factor;
     }
 
-#if 0
-
-    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
-    {
-        OPM_TIMEFUNCTION();
-        // Check if the rates of this well only are single-phase, do nothing
-        // if more than one nonzero rate.
-        auto& ws = well_state.well(this->index_of_well_);
-        // int nonzero_rate_index = -1;
-        // const Scalar floating_point_error_epsilon = 1e-14;
-        // for (int p = 0; p < this->number_of_phases_; ++p) {
-        //     if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
-        //         if (nonzero_rate_index == -1) {
-        //             nonzero_rate_index = p;
-        //         } else {
-        //             // More than one nonzero rate.
-        //             return;
-        //         }
-        //     }
-        // }
-        // Calculate the rates that follow from the current primary variables.
-        std::vector<Scalar> well_q_s;
-        Scalar bhp_tmp = this->isInjector()
-            ? this->well_ecl_.injectionControls(summary_state).bhp_limit
-            : this->well_ecl_.productionControls(summary_state).bhp_limit;
-        computeWellRatesWithBhp(simulator,
-                            bhp_tmp,
-                            well_q_s,
-                            deferred_logger);
-        //                     constcomputeCurrentWellRates(simulator, deferred_logger);
-
-        //if (nonzero_rate_index == -1) {
-            // No nonzero rates.
-            // Use the computed rate directly
-        for (int p = 0; p < this->number_of_phases_; ++p) {
-            ws.surface_rates[p] = well_q_s[this->flowPhaseToModelCompIdx(p)];
-        }
-        // set fractions
-
-        ws.setInitializedFromReservoir(true);
-        //this->updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
-        
-    }
-#endif
 
 
     template <typename TypeTag>
@@ -1659,12 +1608,6 @@ namespace Opm
                         const SummaryState& /* summary_state */,
                         WellState<Scalar>& well_state,
                         DeferredLogger& deferred_logger) const
-    // template <typename TypeTag>
-    // void
-    // WellInterface<TypeTag>::
-    // updateWellStateRates(const Simulator& simulator,
-    //                      WellState<Scalar>& well_state,
-    //                      DeferredLogger& deferred_logger) const
     {
         OPM_TIMEFUNCTION();
         // Check if the rates of this well only are single-phase, do nothing
@@ -1684,8 +1627,10 @@ 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;
+
         if (nonzero_rate_index == -1) {
             // No nonzero rates.
             // Use the computed rate directly
@@ -1694,7 +1639,7 @@ namespace Opm
             }
             return;
         }
-        // set fractions
+
         // 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);