Merge pull request #4686 from GitPaean/support_winjmult_rebase

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.cpp

@@ -1405,4 +1405,28 @@ int BlackoilWellModelGeneric::numLocalNonshutWells() const
     return well_container_generic_.size();
 }
 
+
+void BlackoilWellModelGeneric::initInjMult() {
+    for (auto& well : this->well_container_generic_) {
+        if (well->isInjector() && well->wellEcl().getInjMultMode() != Well::InjMultMode::NONE) {
+            const auto& ws = this->wellState().well(well->indexOfWell());
+            const auto& perf_data = ws.perf_data;
+
+            auto &values = this->prev_inj_multipliers_[well->name()];
+            if (values.empty()) {
+                values.assign(perf_data.size(), 1.0);
+            }
+            well->initInjMult(values);
+        }
+    }
+}
+
+void BlackoilWellModelGeneric::updateInjMult(DeferredLogger& deferred_logger) {
+    for (const auto& well : this->well_container_generic_) {
+        if (well->isInjector() && well->wellEcl().getInjMultMode() != Well::InjMultMode::NONE) {
+            well->updateInjMult(this->prev_inj_multipliers_[well->name()], deferred_logger);
+        }
+    }
+}
+
 }

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -192,6 +192,7 @@ public:
         serializer(closed_this_step_);
         serializer(guideRate_);
         serializer(node_pressures_);
+        serializer(prev_inj_multipliers_);
         serializer(active_wgstate_);
         serializer(last_valid_wgstate_);
         serializer(nupcol_wgstate_);
@@ -208,6 +209,7 @@ public:
                this->local_shut_wells_ == rhs.local_shut_wells_ &&
                this->closed_this_step_ == rhs.closed_this_step_ &&
                this->node_pressures_ == rhs.node_pressures_ &&
+               this->prev_inj_multipliers_ == rhs.prev_inj_multipliers_ &&
                this->active_wgstate_ == rhs.active_wgstate_ &&
                this->last_valid_wgstate_ == rhs.last_valid_wgstate_ &&
                this->nupcol_wgstate_ == rhs.nupcol_wgstate_ &&
@@ -369,6 +371,11 @@ protected:
                               const SummaryConfig& summaryConfig,
                               DeferredLogger& deferred_logger);
 
+    void initInjMult();
+
+
+    void updateInjMult(DeferredLogger& deferred_logger);
+
     // create the well container
     virtual void createWellContainer(const int time_step) = 0;
     virtual void initWellContainer(const int reportStepIdx) = 0;
@@ -428,6 +435,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,11 +310,14 @@ namespace Opm {
             well->setGuideRate(&guideRate_);
         }
 
-        // Close completions due to economical reasons
+        // Close completions due to economic reasons
         for (auto& well : well_container_) {
             well->closeCompletions(wellTestState());
         }
 
+        // we need the inj_multiplier from the previous time step
+        this->initInjMult();
+
         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.
@@ -470,6 +473,10 @@ namespace Opm {
                 well->updateWaterThroughput(dt, this->wellState());
             }
         }
+
+        // at the end of the time step, updating the inj_multiplier saved in WellState for later use
+        this->updateInjMult(local_deferredLogger);
+
         // report well switching
         for (const auto& well : well_container_) {
             well->reportWellSwitching(this->wellState().well(well->indexOfWell()), local_deferredLogger);

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>
@@ -996,7 +997,25 @@ namespace Opm
                               return this->extendEval(value);
                           }
                       };
+
         WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger);
+	
+        if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
+            const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
+            const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index];
+            const int seg = this->segmentNumberToIndex(con.segment());
+            // 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 = segment_pres + 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;
+            }
+        }
     }
 
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -659,8 +659,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,68 @@ 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;
+    // initializing the inj_multipler_ to be 1.0
+    this->inj_multiplier_ = std::vector<double>(max_inj_mult.size(), 1.);
+}
+
+void WellInterfaceGeneric::updateInjMult(std::vector<double>& inj_multipliers, DeferredLogger& deferred_logger) const
+{
+    if (inj_multipliers.size() != this->inj_multiplier_.size()) {
+        OPM_DEFLOG_THROW(std::runtime_error,
+                         fmt::format("We do not support changing connection numbers during simulation with WINJMULT "
+                                     "for well {}", name()),
+                         deferred_logger);
+    }
+
+    inj_multipliers = this->inj_multiplier_;
+}
+
+
+
+double WellInterfaceGeneric::getInjMult(const int perf,
+                                        const double bhp,
+                                        const double perf_pres) const
+{
+    assert(!this->isProducer());
+
+    double multiplier = 1.;
+
+    const auto perf_ecl_index = this->perforationData()[perf].ecl_index;
+    const bool is_wrev = this->well_ecl_.getInjMultMode() == Well::InjMultMode::WREV;
+
+    const bool active_injmult = (is_wrev && this->well_ecl_.aciveWellInjMult()) ||
+                              this->well_ecl_.getConnections()[perf_ecl_index].activeInjMult();
+
+    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;
+
+        const auto frac_press = injmult.fracture_pressure;
+        const auto gradient = injmult.multiplier_gradient;
+        if (pres > frac_press) {
+            multiplier = 1. + (pres - frac_press) * gradient;
+        }
+    }
+
+    // for CIRR mode, if there is no active WINJMULT setup, we will use the previous injection multiplier,
+    // to mimic keeping the existing fracturing open
+    if (this->well_ecl_.getInjMultMode() == Well::InjMultMode::CIRR) {
+        multiplier = std::max(multiplier, this->prev_inj_multiplier_[perf_ecl_index]);
+    }
+
+    this->inj_multiplier_[perf_ecl_index] = multiplier;
+    return multiplier;
+}
+
+
+
+
 bool WellInterfaceGeneric::wellHasTHPConstraints(const SummaryState& summaryState) const
 {
     // only wells under prediction mode can have THP constraint

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -178,6 +178,16 @@ public:
     double wsolvent() const;
     double rsRvInj() const;
 
+    // at the beginning of the time step, we check what inj_multiplier from the previous running
+    void initInjMult(const std::vector<double>& max_inj_mult);
+
+    // update the InjMult information at the end of the time step, so it can be used for later.
+    void updateInjMult(std::vector<double>& inj_multipliers, DeferredLogger& deferred_logger) 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;
 
@@ -352,6 +362,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_;

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;

parallelTests.cmake

@@ -173,6 +173,14 @@ add_test_compare_parallel_simulation(CASENAME actionx_m1
                                      DIR udq_actionx
                                      TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
 
+add_test_compare_parallel_simulation(CASENAME winjmult_msw
+                                     FILENAME WINJMULT_MSW
+                                     SIMULATOR flow
+                                     ABS_TOL ${abs_tol}
+                                     REL_TOL 0.12
+                                     DIR winjmult
+                                     TEST_ARGS --enable-tuning=true)
+
 add_test_compare_parallel_simulation(CASENAME 3_a_mpi_multflt_mod2
                                      FILENAME 3_A_MPI_MULTFLT_SCHED_MODEL2
                                      SIMULATOR flow

regressionTests.cmake

@@ -1176,16 +1176,30 @@ add_test_compareECLFiles(CASENAME 02_wgrupcon
                          ABS_TOL ${abs_tol}
                          REL_TOL ${rel_tol}
                          DIR wgrupcon)
+add_test_compareECLFiles(CASENAME winjmult_stdw
+                         FILENAME WINJMULT_STDW
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR winjmult
+                         TEST_ARGS --enable-tuning=true)
+add_test_compareECLFiles(CASENAME winjmult_msw
+                         FILENAME WINJMULT_MSW
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR winjmult
+                         TEST_ARGS --enable-tuning=true)
 add_test_compareECLFiles(CASENAME 01_vappars
                          FILENAME VAPPARS-01
                          SIMULATOR flow
                          ABS_TOL ${abs_tol}
                          REL_TOL ${rel_tol}
-                         DIR vappars)                         
+                         DIR vappars)
 add_test_compareECLFiles(CASENAME 01_multflt
                          FILENAME MULTFLT-01
                          SIMULATOR flow
                          ABS_TOL ${abs_tol}
                          REL_TOL ${rel_tol}
-                         DIR mult)                         
-                         
+                         DIR mult)
+