Merge pull request #5631 from vkip/winjmult_damp

Add dampening to the injectivity multipliers (WINJMULT)

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -34,6 +34,9 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
 {
     dbhp_max_rel_ = Parameters::Get<Parameters::DbhpMaxRel<Scalar>>();
     dwell_fraction_max_ = Parameters::Get<Parameters::DwellFractionMax<Scalar>>();
+    inj_mult_osc_threshold_ = Parameters::Get<Parameters::InjMultOscThreshold<Scalar>>();
+    inj_mult_damp_mult_ = Parameters::Get<Parameters::InjMultDampMult<Scalar>>();
+    inj_mult_min_damp_factor_ = Parameters::Get<Parameters::InjMultMinDampFactor<Scalar>>();
     max_residual_allowed_ = Parameters::Get<Parameters::MaxResidualAllowed<Scalar>>();
     relaxed_max_pv_fraction_ = Parameters::Get<Parameters::RelaxedMaxPvFraction<Scalar>>();
     tolerance_mb_ = Parameters::Get<Parameters::ToleranceMb<Scalar>>();
@@ -108,6 +111,12 @@ void BlackoilModelParameters<Scalar>::registerParameters()
         ("Maximum relative change of the bottom-hole pressure in a single iteration");
     Parameters::Register<Parameters::DwellFractionMax<Scalar>>
         ("Maximum absolute change of a well's volume fraction in a single iteration");
+    Parameters::Register<Parameters::InjMultOscThreshold<Scalar>>
+        ("Injection multiplier oscillation threshold (used for multiplier dampening)");
+    Parameters::Register<Parameters::InjMultDampMult<Scalar>>
+        ("Injection multiplier dampening factor (dampening multiplied by this each time oscillation is detected)");
+    Parameters::Register<Parameters::InjMultMinDampFactor<Scalar>>
+        ("Minimum injection multiplier dampening factor (maximum dampening level)");
     Parameters::Register<Parameters::MaxResidualAllowed<Scalar>>
         ("Absolute maximum tolerated for residuals without cutting the time step size");
     Parameters::Register<Parameters::RelaxedMaxPvFraction<Scalar>>

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -34,6 +34,15 @@ struct DwellFractionMax { static constexpr Scalar value = 0.2; };
 
 struct EclDeckFileName { static constexpr auto value = ""; };
 
+template<class Scalar>
+struct InjMultOscThreshold { static constexpr Scalar value = 0.1; };
+
+template<class Scalar>
+struct InjMultDampMult { static constexpr Scalar value = 0.9; };
+
+template<class Scalar>
+struct InjMultMinDampFactor { static constexpr Scalar value = 0.05; };
+
 template<class Scalar>
 struct MaxResidualAllowed { static constexpr Scalar value = 1e7; };
 
@@ -156,6 +165,12 @@ public:
     Scalar dbhp_max_rel_;
     /// Max absolute change in well volume fraction in single iteration.
     Scalar dwell_fraction_max_;
+    /// Injectivity multiplier oscillation threshold
+    Scalar inj_mult_osc_threshold_;
+    /// Injectivity multiplier dampening multiplier
+    Scalar inj_mult_damp_mult_;
+    /// Minimum damping factor for injectivity multipliers
+    Scalar inj_mult_min_damp_factor_;
     /// Absolute max limit for residuals.
     Scalar max_residual_allowed_;
     //// Max allowed pore volume faction where CNV is violated. Below the

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -394,6 +394,8 @@ namespace Opm
         // creating a copy of the well itself, to avoid messing up the explicit information
         // during this copy, the only information not copied properly is the well controls
         MultisegmentWell<TypeTag> well_copy(*this);
+        well_copy.resetDampening();
+
         well_copy.debug_cost_counter_ = 0;
 
         // store a copy of the well state, we don't want to update the real well state
@@ -1140,7 +1142,7 @@ namespace Opm
             const Scalar perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
                                                                * this->segments_.perforation_depth_diff(perf);
             const Scalar perf_press = segment_pres + perf_seg_press_diff;
-            const Scalar multiplier = this->getInjMult(perf, segment_pres, perf_press);
+            const Scalar multiplier = this->getInjMult(perf, segment_pres, perf_press, deferred_logger);
             for (std::size_t i = 0; i < mob.size(); ++i) {
                 mob[i] *= multiplier;
             }

opm/simulators/wells/StandardWell_impl.hpp

@@ -689,7 +689,7 @@ namespace Opm
         if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) {
             const Scalar bhp = this->primary_variables_.value(Bhp);
             const Scalar perf_press = bhp +  this->connections_.pressure_diff(perf);
-            const Scalar multiplier = this->getInjMult(perf, bhp, perf_press);
+            const Scalar multiplier = this->getInjMult(perf, bhp, perf_press, deferred_logger);
             for (std::size_t i = 0; i < mob.size(); ++i) {
                 mob[i] *= multiplier;
             }
@@ -1522,6 +1522,7 @@ namespace Opm
         // creating a copy of the well itself, to avoid messing up the explicit information
         // during this copy, the only information not copied properly is the well controls
         StandardWell<TypeTag> well_copy(*this);
+        well_copy.resetDampening();
 
         // iterate to get a more accurate well density
         // create a copy of the well_state to use. If the operability checking is sucessful, we use this one

opm/simulators/wells/WellInterface.hpp

@@ -37,8 +37,6 @@ namespace Opm {
 
 #include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
 
-#include <opm/simulators/flow/BlackoilModelParameters.hpp>
-
 #include <opm/simulators/wells/BlackoilWellModel.hpp>
 #include <opm/simulators/wells/GasLiftGroupInfo.hpp>
 #include <opm/simulators/wells/GasLiftSingleWell.hpp>
@@ -103,7 +101,7 @@ public:
     using WellInterfaceFluidSystem<FluidSystem>::Oil;
     using WellInterfaceFluidSystem<FluidSystem>::Water;
 
-    using ModelParameters = BlackoilModelParameters<Scalar>;
+    using ModelParameters = typename Base::ModelParameters;
 
     static constexpr bool has_solvent = getPropValue<TypeTag, Properties::EnableSolvent>();
     static constexpr bool has_zFraction = getPropValue<TypeTag, Properties::EnableExtbo>();
@@ -385,7 +383,6 @@ public:
                                             const bool fixed_status = false) = 0;
 protected:
     // simulation parameters
-    const ModelParameters& param_;
     std::vector<RateVector> connectionRates_;
     std::vector<Scalar> B_avg_;
     bool changed_to_stopped_this_step_ = false;

opm/simulators/wells/WellInterfaceFluidSystem.cpp

@@ -48,13 +48,14 @@ WellInterfaceFluidSystem<FluidSystem>::
 WellInterfaceFluidSystem(const Well& well,
                          const ParallelWellInfo<Scalar>& parallel_well_info,
                          const int time_step,
+                         const ModelParameters& param,
                          const RateConverterType& rate_converter,
                          const int pvtRegionIdx,
                          const int num_components,
                          const int num_phases,
                          const int index_of_well,
                          const std::vector<PerforationData<Scalar>>& perf_data)
-    : WellInterfaceGeneric<Scalar>(well, parallel_well_info, time_step,
+    : WellInterfaceGeneric<Scalar>(well, parallel_well_info, time_step, param,
                                    pvtRegionIdx, num_components, num_phases,
                                    index_of_well, perf_data)
     , rateConverter_(rate_converter)

opm/simulators/wells/WellInterfaceFluidSystem.hpp

@@ -57,6 +57,7 @@ protected:
 
 public:
     using Scalar = typename FluidSystem::Scalar;
+    using ModelParameters = typename WellInterfaceGeneric<Scalar>::ModelParameters;
 
     int flowPhaseToModelPhaseIdx(const int phaseIdx) const;
 
@@ -73,6 +74,7 @@ protected:
     WellInterfaceFluidSystem(const Well& well,
                              const ParallelWellInfo<Scalar>& parallel_well_info,
                              const int time_step,
+                             const ModelParameters& param,
                              const RateConverterType& rate_converter,
                              const int pvtRegionIdx,
                              const int num_components,

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -57,6 +57,7 @@ WellInterfaceGeneric<Scalar>::
 WellInterfaceGeneric(const Well& well,
                      const ParallelWellInfo<Scalar>& pw_info,
                      const int time_step,
+                     const ModelParameters& param,
                      const int pvtRegionIdx,
                      const int num_components,
                      const int num_phases,
@@ -65,6 +66,7 @@ WellInterfaceGeneric(const Well& well,
       : well_ecl_(well)
       , parallel_well_info_(pw_info)
       , current_step_(time_step)
+      , param_(param)
       , pvtRegionIdx_(pvtRegionIdx)
       , num_components_(num_components)
       , number_of_phases_(num_phases)
@@ -220,8 +222,11 @@ initInjMult(const std::vector<Scalar>& 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
+    const auto nperf = max_inj_mult.size();
-    this->inj_multiplier_ = std::vector<Scalar>(max_inj_mult.size(), 1.);
+    // initializing the (report?) step specific multipliers and damping factors to be 1.0
+    this->inj_multiplier_ = std::vector<Scalar>(nperf, 1.);
+    this->inj_multiplier_previter_ = std::vector<Scalar>(nperf, 1.);
+    this->inj_multiplier_damp_factor_ = std::vector<Scalar>(nperf, 1.);
 }
 
 template<class Scalar>
@@ -243,7 +248,8 @@ template<class Scalar>
 Scalar WellInterfaceGeneric<Scalar>::
 getInjMult(const int perf,
            const Scalar bhp,
-           const Scalar perf_pres) const
+           const Scalar perf_pres,
+           DeferredLogger& dlogger) const
 {
     assert(!this->isProducer());
 
@@ -265,8 +271,29 @@ getInjMult(const int perf,
         if (pres > frac_press) {
             multiplier = 1. + (pres - frac_press) * gradient;
         }
-    }
+        const Scalar osc_threshold = param_.inj_mult_osc_threshold_;
+        const Scalar prev_multiplier = this->inj_multiplier_[perf_ecl_index] > 0 ? this->inj_multiplier_[perf_ecl_index] : 1.0;
+        if (std::abs(multiplier - prev_multiplier) > osc_threshold) {
+            const Scalar prev2_multiplier = this->inj_multiplier_previter_[perf_ecl_index] > 0 ? this->inj_multiplier_previter_[perf_ecl_index] : 1.0;
+            const bool oscillating = (multiplier - prev_multiplier) * (prev_multiplier - prev2_multiplier) < 0;
 
+            const Scalar min_damp_factor = this->param_.inj_mult_min_damp_factor_;
+            Scalar damp_factor = this->inj_multiplier_damp_factor_[perf_ecl_index];
+            if (oscillating) {
+                if (damp_factor > min_damp_factor) {
+                    const Scalar damp_multiplier = this->param_.inj_mult_damp_mult_;
+                    damp_factor = std::max(min_damp_factor, damp_factor * damp_multiplier);
+                    this->inj_multiplier_damp_factor_[perf_ecl_index] = damp_factor;
+                } else {
+                    const auto msg = fmt::format("Well {}, perforation {}: Injectivity multiplier dampening factor reached minimum (= {})",
+                                                 this->name(),
+                                                 perf_ecl_index,
+                                                 min_damp_factor);
+                    dlogger.warning(msg);
+                }
+            }
+            multiplier = multiplier * damp_factor + (1.0 - damp_factor) * prev_multiplier;
+        }
         // 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) {
@@ -274,6 +301,9 @@ getInjMult(const int perf,
         }
 
         this->inj_multiplier_[perf_ecl_index] = multiplier;
+        this->inj_multiplier_previter_[perf_ecl_index] = prev_multiplier;
+    }
+
     return multiplier;
 }
 

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -25,6 +25,7 @@
 #define OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED
 
 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
+#include <opm/simulators/flow/BlackoilModelParameters.hpp>
 
 #include <map>
 #include <optional>
@@ -50,9 +51,12 @@ class Schedule;
 template<class Scalar>
 class WellInterfaceGeneric {
 public:
+    using ModelParameters = BlackoilModelParameters<Scalar>;
+
     WellInterfaceGeneric(const Well& well,
                          const ParallelWellInfo<Scalar>& parallel_well_info,
                          const int time_step,
+                         const ModelParameters& param,
                          const int pvtRegionIdx,
                          const int num_components,
                          const int num_phases,
@@ -156,7 +160,7 @@ public:
 
     // Note:: for multisegment wells, bhp is actually segment pressure in practice based on observation
     // it might change in the future
-    Scalar getInjMult(const int perf, const Scalar bhp, const Scalar perf_pres) const;
+    Scalar getInjMult(const int perf, const Scalar bhp, const Scalar perf_pres, DeferredLogger& dlogger) const;
 
     // whether a well is specified with a non-zero and valid VFP table number
     bool isVFPActive(DeferredLogger& deferred_logger) const;
@@ -217,6 +221,10 @@ protected:
                                          Well::InjectionControls& inj_controls,
                                          Well::ProductionControls& prod_controls) const;
 
+    void resetDampening() {
+        std::fill(this->inj_multiplier_damp_factor_.begin(), this->inj_multiplier_damp_factor_.end(), 1.0);
+    }
+
     // definition of the struct OperabilityStatus
     struct OperabilityStatus
     {
@@ -274,6 +282,7 @@ protected:
 
     const ParallelWellInfo<Scalar>& parallel_well_info_;
     const int current_step_;
+    const ModelParameters& param_;
 
     // The pvt region of the well. We assume
     // We assume a well to not penetrate more than one pvt region.
@@ -355,6 +364,11 @@ protected:
     // which intends to keep the fracturing open
     std::vector<Scalar> prev_inj_multiplier_;
 
+    // WINJMULT multipliers for previous iteration (used for oscillation detection)
+    mutable std::vector<Scalar> inj_multiplier_previter_;
+    // WINJMULT dampening factors (used in case of oscillations)
+    mutable std::vector<Scalar> inj_multiplier_damp_factor_;
+
     // the multiplier due to injection filtration cake
     std::vector<Scalar> inj_fc_multiplier_;
 

opm/simulators/wells/WellInterfaceIndices.cpp

@@ -38,6 +38,7 @@ WellInterfaceIndices<FluidSystem,Indices>::
 WellInterfaceIndices(const Well& well,
                      const ParallelWellInfo<Scalar>& parallel_well_info,
                      const int time_step,
+                     const ModelParameters& param,
                      const typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType& rate_converter,
                      const int pvtRegionIdx,
                      const int num_components,
@@ -47,6 +48,7 @@ WellInterfaceIndices(const Well& well,
     : WellInterfaceFluidSystem<FluidSystem>(well,
                                             parallel_well_info,
                                             time_step,
+                                            param,
                                             rate_converter,
                                             pvtRegionIdx,
                                             num_components,

opm/simulators/wells/WellInterfaceIndices.hpp

@@ -38,6 +38,7 @@ public:
     using WellInterfaceFluidSystem<FluidSystem>::Water;
     using Scalar = typename FluidSystem::Scalar;
     using Eval = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq>;
+    using ModelParameters = typename WellInterfaceFluidSystem<FluidSystem>::ModelParameters;
 
     int flowPhaseToModelCompIdx(const int phaseIdx) const;
     int modelCompIdxToFlowCompIdx(const unsigned compIdx) const;
@@ -58,6 +59,7 @@ protected:
     WellInterfaceIndices(const Well& well,
                          const ParallelWellInfo<Scalar>& parallel_well_info,
                          const int time_step,
+                         const ModelParameters& param,
                          const typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType& rate_converter,
                          const int pvtRegionIdx,
                          const int num_components,

opm/simulators/wells/WellInterface_impl.hpp

@@ -65,13 +65,13 @@ namespace Opm
       : WellInterfaceIndices<FluidSystem,Indices>(well,
                                                   pw_info,
                                                   time_step,
+                                                  param,
                                                   rate_converter,
                                                   pvtRegionIdx,
                                                   num_components,
                                                   num_phases,
                                                   index_of_well,
                                                   perf_data)
-      , param_(param)
     {
         connectionRates_.resize(this->number_of_perforations_);
 
@@ -212,17 +212,17 @@ namespace Opm
         } else {
             from = WellProducerCMode2String(ws.production_cmode);
         }
-        bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
+        bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
 
         if (oscillating) {
             // only output frist time
-            bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == param_.max_number_of_well_switches_;
+            bool output = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) == this->param_.max_number_of_well_switches_;
             if (output) {
                 std::ostringstream ss;
                 ss << "    The control mode for well " << this->name()
                    << " is oscillating\n"
                    << "    We don't allow for more than "
-                   << param_.max_number_of_well_switches_
+                   << this->param_.max_number_of_well_switches_
                    << " switches. The control is kept at " << from;
                 deferred_logger.info(ss.str());
                 // add one more to avoid outputting the same info again
@@ -287,7 +287,7 @@ namespace Opm
         } else {
             from = WellProducerCMode2String(ws.production_cmode);
         }
-        const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= param_.max_number_of_well_switches_;
+        const bool oscillating = std::count(this->well_control_log_.begin(), this->well_control_log_.end(), from) >= this->param_.max_number_of_well_switches_;
 
         if (oscillating || this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger) || !(this->well_ecl_.getStatus() == WellStatus::OPEN)) {
            return false;
@@ -309,7 +309,7 @@ namespace Opm
                                                                       prod_controls.hasControl(Well::ProducerCMode::GRUP);
 
                     changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls);
-                    if (hasGroupControl && param_.check_group_constraints_inner_well_iterations_) {
+                    if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) {
                         changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state,deferred_logger);
                     }
 
@@ -711,7 +711,7 @@ namespace Opm
             deferred_logger.debug("WellTest: Well equation for well " + this->name() +  " converged");
             return true;
         }
-        const int max_iter = param_.max_welleq_iter_;
+        const int max_iter = this->param_.max_welleq_iter_;
         deferred_logger.debug("WellTest: Well equation for well " + this->name() + " failed converging in "
                               + std::to_string(max_iter) + " iterations");
         well_state = well_state0;
@@ -767,7 +767,7 @@ namespace Opm
         }
 
         if (!converged) {
-            const int max_iter = param_.max_welleq_iter_;
+            const int max_iter = this->param_.max_welleq_iter_;
             deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "
                                   + std::to_string(max_iter) + " iterations");
             well_state = well_state0;
@@ -821,18 +821,18 @@ namespace Opm
     {
         const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
 
-        if (param_.check_well_operability_iter_)
+        if (this->param_.check_well_operability_iter_)
             checkWellOperability(simulator, well_state, deferred_logger);
 
         // only use inner well iterations for the first newton iterations.
         const int iteration_idx = simulator.model().newtonMethod().numIterations();
-        if (iteration_idx < param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
+        if (iteration_idx < this->param_.max_niter_inner_well_iter_ || this->well_ecl_.isMultiSegment()) {
             this->operability_status_.solvable = true;
             bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger);
 
             // unsolvable wells are treated as not operable and will not be solved for in this iteration.
             if (!converged) {
-                if (param_.shut_unsolvable_wells_)
+                if (this->param_.shut_unsolvable_wells_)
                     this->operability_status_.solvable = false;
             }
         }
@@ -912,7 +912,7 @@ namespace Opm
                          DeferredLogger& deferred_logger)
     {
 
-        if (!param_.check_well_operability_) {
+        if (!this->param_.check_well_operability_) {
             return;
         }