store the inj_multiper in PerfData in WellState

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -428,9 +428,6 @@ protected:
     std::unique_ptr<VFPProperties> vfp_properties_{};
     std::map<std::string, double> node_pressures_; // Storing network pressures for output.
 
-    // previous injection multiplier, it is used in the injection multiplier calculation for WINJMULT keyword
-    std::unordered_map<std::string, std::vector<double>> prev_inj_multipliers_;
-
     /*
       The various wellState members should be accessed and modified
       through the accessor functions wellState(), prevWellState(),

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -310,15 +310,12 @@ namespace Opm {
             well->setGuideRate(&guideRate_);
         }
 
-        // initialize the WINJMULT related, we need to record the highest injecting multiplier historically
+        // we need the inj_multiplier from the previous time step
-        // to support the CIRR mode of  WINJMULT keyword
         for (auto& well : well_container_) {
             if (well->isInjector()) {
-                auto& ws = this->wellState().well(well->indexOfWell());
+                const auto& ws = this->wellState().well(well->indexOfWell());
-                if ((this->prev_inj_multipliers_.count(well->name())) == 0 ) {
+                const auto& perf_data = ws.perf_data;
-                    this->prev_inj_multipliers_[well->name()] = std::vector<double>(ws.perf_data.size(), 1.0);
+                well->initInjMult(perf_data.inj_multiplier);
-                }
-                well->initInjMult(this->prev_inj_multipliers_.at(well->name()));
             }
         }
 
@@ -481,8 +478,10 @@ namespace Opm {
             if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
                 well->updateWaterThroughput(dt, this->wellState());
             }
+            // at the end of the time step, updating the inj_multiplier saved in WellState for later use
             if (well->isInjector()) {
-                well->updateMaxInjMult(this->prev_inj_multipliers_[well->name()]);
+                auto& perf_data = this->wellState().well(well->indexOfWell()).perf_data;
+                well->updateInjMult(perf_data.inj_multiplier);
             }
         }
         // report well switching

opm/simulators/wells/PerfData.cpp

@@ -43,6 +43,7 @@ PerfData::PerfData(std::size_t num_perf, double pressure_first_connection_, bool
     , ecl_index(num_perf)
 {
     if (injector) {
+        this->inj_multiplier.resize(num_perf, 1.0);
         this->water_throughput.resize(num_perf);
         this->skin_pressure.resize(num_perf);
         this->water_velocity.resize(num_perf);
@@ -94,6 +95,7 @@ bool PerfData::try_assign(const PerfData& other) {
     this->solvent_rates = other.solvent_rates;
     this->polymer_rates = other.polymer_rates;
     this->brine_rates = other.brine_rates;
+    this->inj_multiplier = other.inj_multiplier;
     this->water_throughput = other.water_throughput;
     this->skin_pressure = other.skin_pressure;
     this->water_velocity = other.water_velocity;
@@ -117,6 +119,7 @@ bool PerfData::operator==(const PerfData& rhs) const
            this->connection_transmissibility_factor == rhs.connection_transmissibility_factor &&
            this->satnum_id == rhs.satnum_id &&
            this->ecl_index == rhs.ecl_index &&
+           this->inj_multiplier == rhs.inj_multiplier &&
            this->water_throughput == rhs.water_throughput &&
            this->skin_pressure == rhs.skin_pressure &&
            this->water_velocity == rhs.water_velocity;

opm/simulators/wells/PerfData.hpp

@@ -57,6 +57,7 @@ public:
         serializer(connection_transmissibility_factor);
         serializer(satnum_id);
         serializer(ecl_index);
+        serializer(inj_multiplier);
         serializer(water_throughput);
         serializer(skin_pressure);
         serializer(water_velocity);
@@ -78,6 +79,9 @@ public:
     std::vector<double> connection_transmissibility_factor;
     std::vector<int> satnum_id;
     std::vector<std::size_t> ecl_index;
+    // the injection multiplier due to WINJMULT keyword
+    // it only applies to injectors
+    std::vector<double> inj_multiplier;
 
     // The water_throughput, skin_pressure and water_velocity variables are only
     // used for injectors to check the injectivity.

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -195,15 +195,15 @@ void WellInterfaceGeneric::initInjMult(const std::vector<double>& max_inj_mult)
     // prev_inj_multiplier_ will stay unchanged during the time step
     // while inj_multiplier_ might be updated during the time step
     this->prev_inj_multiplier_ = max_inj_mult;
-    // reset the inj_multipler_ to be 1.0
+    // initializing the inj_multipler_ to be 1.0
     this->inj_multiplier_ = std::vector<double>(max_inj_mult.size(), 1.);
 }
 
-void WellInterfaceGeneric::updateMaxInjMult(std::vector<double>& max_multipliers) const
+void WellInterfaceGeneric::updateInjMult(std::vector<double>& inj_multipliers) const
 {
-    assert(max_multipliers.size() == this->inj_multiplier_.size());
+    assert(inj_multipliers.size() == this->inj_multiplier_.size());
 
-    max_multipliers = this->inj_multiplier_;
+    inj_multipliers = this->inj_multiplier_;
 }
 
 
@@ -214,16 +214,19 @@ double WellInterfaceGeneric::getInjMult(const int perf,
 {
     assert(!this->isProducer());
 
+    double multipler = 1.;
+
     const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
     const bool is_wrev = this->well_ecl_.getInjMultMode() == Well::InjMultMode::WREV;
 
-    const auto& injmult = is_wrev ?  this->well_ecl_.getWellInjMult() :
+    const bool active_injmult = (is_wrev && this->well_ecl_.aciveWellInjMult()) ||
-                                                  this->well_ecl_.getConnections()[perf_ecl_index].injmult();
+                              this->well_ecl_.getConnections()[perf_ecl_index].activeInjMult();
-    const double pres = is_wrev ? bhp : perf_pres;
 
+    if (active_injmult) {
+        const auto& injmult= is_wrev ? this->well_ecl_.getWellInjMult()
+                                                  : this->well_ecl_.getConnections()[perf_ecl_index].injmult();
+        const double pres = is_wrev ? bhp : perf_pres;
 
-    double multipler = 1.;
-    if (injmult.active()) {
         const auto frac_press = injmult.fracture_pressure;
         const auto gradient = injmult.multiplier_gradient;
         if (pres > frac_press) {
@@ -231,6 +234,8 @@ double WellInterfaceGeneric::getInjMult(const int perf,
         }
     }
 
+    // for CIRR mode, if there is not activel WINJMULT setup, we will use the previous injection multiplier,
+    // to mimic keeping the existing fracturing open
     if (this->well_ecl_.getInjMultMode() == Well::InjMultMode::CIRR) {
         multipler = std::max(multipler, this->prev_inj_multiplier_[perf_ecl_index]);
     }

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -178,16 +178,11 @@ public:
     double wsolvent() const;
     double rsRvInj() const;
 
-    // when creating the well, setting the maximum injection multiplier
+    // at the beginning of the time step, we check what inj_multiplier from the previous running
-    // it can be used in the multiplier calculation with keyword WINJMULT.
-    // and also reset the inj_multiplier_ to be 1 to start.
-    // the reason we need to have two sets of values is because for CIRR mode,
-    // we will only update the maximum multiplier achieved so far after the solution gets converged,
-    // which gives the best results during testing.
     void initInjMult(const std::vector<double>& max_inj_mult);
 
-    // update the InjMult information in the BlackoilWellModel at the end of the time step
+    // update the InjMult information at the end of the time step, so it can be used for later.
-    void updateMaxInjMult(std::vector<double>& max_multipliers) const;
+    void updateInjMult(std::vector<double>& inj_multipliers) const;
 
     // Note:: for multisegment wells, bhp is actually segment pressure in practice based on observation
     // it might change in the future

opm/simulators/wells/WellState.cpp

@@ -378,6 +378,7 @@ void WellState::init(const std::vector<double>& cellPressures,
                 // perforations is equal, otherwise initialize
                 // perfphaserates to well rates divided by the
                 // number of perforations.
+                // TODO: we might still need the values from the prev_well if the connection structure changes
                 if (global_num_perf_same)
                 {
                     auto& perf_data = new_well.perf_data;