Merge 3e3b72d037 into 4b688e5945

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -436,10 +436,10 @@ 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);
+                //const bool zero_target = well->stoppedOrZeroRateTarget(simulator_, this->wellState(), local_deferredLogger);
-                if (well->isProducer() && !zero_target) {
+                //if (well->isProducer()){// && !zero_target) {
-                    well->updateWellStateRates(simulator_, this->wellState(), local_deferredLogger);
+                    well->initializeWellState(simulator_, this->groupState(), this->wellState(), local_deferredLogger);
-                }
+                //}
             }
         }
 
@@ -559,7 +559,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->initializeWellState(simulator_, this->groupState(), this->wellState(), deferred_logger);
             }
             if (well->isVFPActive(deferred_logger)) {
                 well->setPrevSurfaceRates(this->wellState(), this->prevWellState());

opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp

@@ -111,6 +111,18 @@ 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.stw_primaryvar.size() == value_.size();
+        if(prim_set){
+        //if (std::abs(total_seg_rate) > 0.) {
+            if (has_wfrac_variable) {
+                value_[seg][WFrac] = ws.multiseg_primaryvar[seg][WFrac];
+            }
+            if (has_gfrac_variable) {
+                value_[seg][GFrac] = ws.multiseg_primaryvar[seg][GFrac];
+            }
+        } else {
         if (std::abs(total_seg_rate) > 0.) {
             if (has_wfrac_variable) {
                 const int water_pos = pu.phase_pos[Water];
@@ -120,6 +132,7 @@ update(const WellState<Scalar>& well_state,
                 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; 
             }
+            // what about water and gas injection?
          } else { // total_seg_rate == 0
             if (well_.isInjector()) {
                 // only single phase injection handled
@@ -152,6 +165,7 @@ update(const WellState<Scalar>& well_state,
             }
         }
         }
+    }
     setEvaluationsFromValues();
 }
 
@@ -230,6 +244,9 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
     const int oil_pos = pu.phase_pos[Oil];
 
     auto& ws = well_state.well(well_.indexOfWell());
+    if constexpr ( numWellEq == 4) {
+        ws.multiseg_primaryvar = value_;
+    }
     auto& segments = ws.segments;
     auto& segment_rates = segments.rates;
     auto& disgas = segments.dissolved_gas_rate;

opm/simulators/wells/SingleWellState.cpp

@@ -51,6 +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_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;

opm/simulators/wells/SingleWellState.hpp

@@ -76,6 +76,9 @@ public:
         serializer(production_cmode);
         serializer(filtrate_conc);
         serializer(perf_data);
+        serializer(stw_primaryvar);
+        serializer(multiseg_primaryvar);
+        serializer(initialized_from_reservoir_);
     }
 
     bool operator==(const SingleWellState&) const;
@@ -142,8 +145,13 @@ 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;
+
     Scalar sum_connection_rates(const std::vector<Scalar>& connection_rates) const;
 };
 

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -162,7 +162,22 @@ 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){   
+        if constexpr (has_wfrac_variable) {
+             value_[WFrac] = ws.stw_primaryvar[WFrac];
+        }
+        if constexpr (has_gfrac_variable) {
+             value_[GFrac] = ws.stw_primaryvar[GFrac];
+        }
+        if constexpr (Indices::enableSolvent) {
+            value_[SFrac] = ws.stw_primaryvar[SFrac];
+         }
+     }else{
+     assert(prim_set==false);   
      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 +190,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
@@ -220,8 +234,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();
 }
@@ -334,7 +348,6 @@ 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;
@@ -399,6 +412,7 @@ copyToWellState(WellState<Scalar>& well_state,
     }
     
     auto& ws = well_state.well(well_.indexOfWell());
+    ws.stw_primaryvar = value_;
     ws.bhp = value_[Bhp];
 
     // calculate the phase rates based on the primary variables

opm/simulators/wells/WellInterface.hpp

@@ -340,10 +340,15 @@ 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,
+    void initializeWellState(const Simulator& simulator,
+                             const GroupState<Scalar>& group_state,
                              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

@@ -1621,7 +1621,8 @@ namespace Opm
     template <typename TypeTag>
     void
     WellInterface<TypeTag>::
-    updateWellStateRates(const Simulator& simulator,
+    initializeWellState(const Simulator& simulator,
+                        const GroupState<Scalar>& group_state,
                         WellState<Scalar>& well_state,
                         DeferredLogger& deferred_logger) const
     {
@@ -1629,44 +1630,90 @@ namespace Opm
         // 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;
+        // int nonzero_rate_index = -1;
-        const Scalar floating_point_error_epsilon = 1e-14;
+        // const Scalar floating_point_error_epsilon = 1e-14;
-        for (int p = 0; p < this->number_of_phases_; ++p) {
+        // for (int p = 0; p < this->number_of_phases_; ++p) {
-            if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
+        //     if (std::abs(ws.surface_rates[p]) > floating_point_error_epsilon) {
-                if (nonzero_rate_index == -1) {
+        //         if (nonzero_rate_index == -1) {
-                    nonzero_rate_index = p;
+        //             nonzero_rate_index = p;
-                } else {
+        //         } else {
-                    // More than one nonzero rate.
+        //             // More than one nonzero rate.
-                    return;
+        //             return;
-                }
+        //         }
-            }
+        //     }
-        }
+        // }
-
         // Calculate the rates that follow from the current primary variables.
-        std::vector<Scalar> well_q_s = computeCurrentWellRates(simulator, deferred_logger);
+        std::vector<Scalar> well_q_s;
+        Scalar bhp_tmp = 0.0;
+        computeWellRatesWithBhp(simulator,
+                            bhp_tmp,
+                            well_q_s,
+                            deferred_logger);
+        //                     constcomputeCurrentWellRates(simulator, deferred_logger);
 
-        if (nonzero_rate_index == -1) {
+        //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)];
         }
-            return;
+        // set fractions
+
+        ws.setInitializedFromReservoir(true);
+        this->updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger);
+        
     }
 
-        // 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];
+    // template <typename TypeTag>
-        const int comp_idx_nz = this->flowPhaseToModelCompIdx(nonzero_rate_index);
+    // void
-        if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) {
+    // WellInterface<TypeTag>::
-            for (int p = 0; p < this->number_of_phases_; ++p) {
+    // updateWellStateRates(const Simulator& simulator,
-                if (p != nonzero_rate_index) {
+    //                      WellState<Scalar>& well_state,
-                    const int comp_idx = this->flowPhaseToModelCompIdx(p);
+    //                      DeferredLogger& deferred_logger) const
-                    Scalar& rate = ws.surface_rates[p];
+    // {
-                    rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
+    //     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 = computeCurrentWellRates(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)];
+    //         }
+    //         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);
+    //     if (std::abs(well_q_s[comp_idx_nz]) > 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]);
+    //             }
+    //         }
+    //     }
+    // }
 
     template <typename TypeTag>
     std::vector<typename WellInterface<TypeTag>::Scalar>