support WINJMULT

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -711,7 +711,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"WH3NUM", {true, std::nullopt}},
         {"WHEDREFD", {true, std::nullopt}},
         {"WHTEMP", {true, std::nullopt}},
-        {"WINJMULT", {true, std::nullopt}},
         {"WLIMTOL", {true, std::nullopt}},
         {"WLIFT", {true, std::nullopt}},
         {"WLISTARG", {true, std::nullopt}},

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -428,6 +428,9 @@ 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,7 +310,19 @@ namespace Opm {
             well->setGuideRate(&guideRate_);
         }
 
-        // Close completions due to economical reasons
+        // initialize the WINJMULT related, we need to record the highest injecting multiplier historically
+        // to support the CIRR mode of  WINJMULT keyword
+        for (auto& well : well_container_) {
+            if (well->isInjector()) {
+                auto& ws = this->wellState().well(well->indexOfWell());
+                if ((this->prev_inj_multipliers_.count(well->name())) == 0 ) {
+                    this->prev_inj_multipliers_[well->name()] = std::vector<double>(ws.perf_data.size(), 1.0);
+                }
+                well->initInjMult(this->prev_inj_multipliers_.at(well->name()));
+            }
+        }
+
+        // Close completions due to economic reasons
         for (auto& well : well_container_) {
             well->closeCompletions(wellTestState());
         }
@@ -469,6 +481,9 @@ namespace Opm {
             if (getPropValue<TypeTag, Properties::EnablePolymerMW>() && well->isInjector()) {
                 well->updateWaterThroughput(dt, this->wellState());
             }
+            if (well->isInjector()) {
+                well->updateMaxInjMult(this->prev_inj_multipliers_[well->name()]);
+            }
         }
         // report well switching
         for (const auto& well : well_container_) {

opm/simulators/wells/MultisegmentWell.hpp

@@ -220,6 +220,7 @@ namespace Opm
         // get the mobility for specific perforation
         template<class Value>
         void getMobility(const Simulator& ebosSimulator,
+                         const int seg,
                          const int perf,
                          std::vector<Value>& mob,
                          DeferredLogger& deferred_logger) const;

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -22,6 +22,7 @@
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/input/eclipse/Schedule/MSW/Valve.hpp>
+#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 #include <opm/input/eclipse/Units/Units.hpp>
 
 #include <opm/material/densead/EvaluationFormat.hpp>
@@ -340,7 +341,7 @@ namespace Opm
                 const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
                 // flux for each perforation
                 std::vector<Scalar> mob(this->num_components_, 0.);
-                getMobility(ebosSimulator, perf, mob, deferred_logger);
+                getMobility(ebosSimulator, seg, perf, mob, deferred_logger);
                 double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
                 const double Tw = this->well_index_[perf] * trans_mult;
 
@@ -669,7 +670,10 @@ namespace Opm
             };
 
             std::vector<Scalar> mob(this->num_components_, 0.0);
-            getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
+            const WellConnections& connections = this->wellEcl().getConnections();
+            const Connection& connection = connections.get(static_cast<int>(subsetPerfID));
+            const int seg = this->segmentNumberToIndex(connection.segment());
+            getMobility(ebosSimulator, seg, static_cast<int>(subsetPerfID), mob, deferred_logger);
 
             const auto& fs = fluidState(subsetPerfID);
             setToZero(connPI);
@@ -982,6 +986,7 @@ namespace Opm
     void
     MultisegmentWell<TypeTag>::
     getMobility(const Simulator& ebosSimulator,
+                const int seg,
                 const int perf,
                 std::vector<Value>& mob,
                 DeferredLogger& deferred_logger) const
@@ -996,6 +1001,20 @@ namespace Opm
                           }
                       };
         WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger);
+	
+        if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
+            // from the reference results, it looks like MSW uses segment pressure instead of BHP here
+            // Note: this is against the documented definition.
+            // we can change this depending on what we want
+            const double segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
+            const double perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
+                                                               * this->segments_.perforation_depth_diff(perf);
+            const double perf_press = this->primary_variables_.getSegmentPressure(seg).value() + perf_seg_press_diff;
+            const double multiplier = this->getInjMult(perf, segment_pres, perf_press);
+            for (size_t i = 0; i < mob.size(); ++i) {
+                mob[i] *= multiplier;
+            }
+        }
     }
 
 
@@ -1095,7 +1114,7 @@ namespace Opm
                 std::vector<Scalar> mob(this->num_components_, 0.0);
 
                 // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
-                getMobility(ebos_simulator, perf, mob, deferred_logger);
+                getMobility(ebos_simulator, seg, perf, mob, deferred_logger);
 
                 const int cell_idx = this->well_cells_[perf];
                 const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
@@ -1453,7 +1472,7 @@ namespace Opm
                 const int cell_idx = this->well_cells_[perf];
                 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
                 std::vector<EvalWell> mob(this->num_components_, 0.0);
-                getMobility(ebosSimulator, perf, mob, deferred_logger);
+                getMobility(ebosSimulator, seg, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
                 const double Tw = this->well_index_[perf] * trans_mult;
                 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
@@ -1765,7 +1784,7 @@ namespace Opm
                 const int cell_idx = this->well_cells_[perf];
                 const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
                 std::vector<Scalar> mob(this->num_components_, 0.0);
-                getMobility(ebosSimulator, perf, mob, deferred_logger);
+                getMobility(ebosSimulator, seg, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
                 const double Tw = this->well_index_[perf] * trans_mult;
                 std::vector<Scalar> cq_s(this->num_components_, 0.0);

opm/simulators/wells/StandardWell_impl.hpp

@@ -610,8 +610,17 @@ namespace Opm
                 }
             }
         }
-    }
 
+        // if the injecting well has WINJMULT setup, we update the mobility accordingly
+        if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
+            const double bhp = this->primary_variables_.value(Bhp);
+            const double perf_press = bhp +  this->connections_.pressure_diff(perf);
+            const double multiplier = this->getInjMult(perf, bhp, perf_press);
+            for (size_t i = 0; i < mob.size(); ++i) {
+                mob[i] *= multiplier;
+            }
+        }
+    }
 
 
     template<typename TypeTag>

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -190,6 +190,58 @@ double WellInterfaceGeneric::rsRvInj() const
     return well_ecl_.getInjectionProperties().rsRvInj;
 }
 
+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
+    this->inj_multiplier_ = std::vector<double>(max_inj_mult.size(), 1.);
+}
+
+void WellInterfaceGeneric::updateMaxInjMult(std::vector<double>& max_multipliers) const
+{
+    assert(max_multipliers.size() == this->inj_multiplier_.size());
+
+    max_multipliers = this->inj_multiplier_;
+}
+
+
+
+double WellInterfaceGeneric::getInjMult(const int perf,
+                                        const double bhp,
+                                        const double perf_pres) const
+{
+    assert(!this->isProducer());
+
+    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() :
+                                                  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) {
+            multipler = 1. + (pres - frac_press) * gradient;
+        }
+    }
+
+    if (this->well_ecl_.getInjMultMode() == Well::InjMultMode::CIRR) {
+        multipler = std::max(multipler, this->prev_inj_multiplier_[perf_ecl_index]);
+    }
+
+    this->inj_multiplier_[perf_ecl_index] = multipler;
+    return multipler;
+}
+
+
+
+
 bool WellInterfaceGeneric::wellHasTHPConstraints(const SummaryState& summaryState) const
 {
     // only wells under prediction mode can have THP constraint

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -178,6 +178,21 @@ public:
     double wsolvent() const;
     double rsRvInj() const;
 
+    // when creating the well, setting the maximum injection multiplier
+    // 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
+    void updateMaxInjMult(std::vector<double>& max_multipliers) const;
+
+    // Note:: for multisegment wells, bhp is actually segment pressure in practice based on observation
+    // it might change in the future
+    double getInjMult(const int perf, const double bhp, const double perf_pres) const;
+
     // whether a well is specified with a non-zero and valid VFP table number
     bool isVFPActive(DeferredLogger& deferred_logger) const;
 
@@ -348,6 +363,13 @@ protected:
     double wsolvent_;
     std::optional<double> dynamic_thp_limit_;
 
+    // recording the multiplier calculate from the keyword WINJMULT during the time step
+    mutable std::vector<double> inj_multiplier_;
+
+    // the injection multiplier from the previous running, it is mostly used for CIRR mode
+    // which intends to keep the fracturing open
+    std::vector<double> prev_inj_multiplier_;
+
     double well_efficiency_factor_;
     const VFPProperties* vfp_properties_;
     const GuideRate* guide_rate_;