diff --git a/opm/simulators/flow/BlackoilModelParameters.cpp b/opm/simulators/flow/BlackoilModelParameters.cpp index 9165c813e..95b26ff0b 100644 --- a/opm/simulators/flow/BlackoilModelParameters.cpp +++ b/opm/simulators/flow/BlackoilModelParameters.cpp @@ -34,6 +34,9 @@ BlackoilModelParameters::BlackoilModelParameters() { dbhp_max_rel_ = Parameters::Get>(); dwell_fraction_max_ = Parameters::Get>(); + inj_mult_osc_threshold_ = Parameters::Get>(); + inj_mult_damp_mult_ = Parameters::Get>(); + inj_mult_min_damp_factor_ = Parameters::Get>(); max_residual_allowed_ = Parameters::Get>(); relaxed_max_pv_fraction_ = Parameters::Get>(); tolerance_mb_ = Parameters::Get>(); @@ -108,6 +111,12 @@ void BlackoilModelParameters::registerParameters() ("Maximum relative change of the bottom-hole pressure in a single iteration"); Parameters::Register> ("Maximum absolute change of a well's volume fraction in a single iteration"); + Parameters::Register> + ("Injection multiplier oscillation threshold (used for multiplier dampening)"); + Parameters::Register> + ("Injection multiplier dampening factor (dampening multiplied by this each time oscillation is detected)"); + Parameters::Register> + ("Minimum injection multiplier dampening factor (maximum dampening level)"); Parameters::Register> ("Absolute maximum tolerated for residuals without cutting the time step size"); Parameters::Register> diff --git a/opm/simulators/flow/BlackoilModelParameters.hpp b/opm/simulators/flow/BlackoilModelParameters.hpp index 524d276db..d654a6fa5 100644 --- a/opm/simulators/flow/BlackoilModelParameters.hpp +++ b/opm/simulators/flow/BlackoilModelParameters.hpp @@ -34,6 +34,15 @@ struct DwellFractionMax { static constexpr Scalar value = 0.2; }; struct EclDeckFileName { static constexpr auto value = ""; }; +template +struct InjMultOscThreshold { static constexpr Scalar value = 0.1; }; + +template +struct InjMultDampMult { static constexpr Scalar value = 0.9; }; + +template +struct InjMultMinDampFactor { static constexpr Scalar value = 0.05; }; + template 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 diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 68e771e45..8ab93a5e1 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1140,7 +1140,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; } diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 95404ddd4..72c047b4b 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/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; } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 32f69e302..7257f3846 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -37,8 +37,6 @@ namespace Opm { #include -#include - #include #include #include @@ -103,7 +101,7 @@ public: using WellInterfaceFluidSystem::Oil; using WellInterfaceFluidSystem::Water; - using ModelParameters = BlackoilModelParameters; + using ModelParameters = typename Base::ModelParameters; static constexpr bool has_solvent = getPropValue(); static constexpr bool has_zFraction = getPropValue(); @@ -385,7 +383,6 @@ public: const bool fixed_status = false) = 0; protected: // simulation parameters - const ModelParameters& param_; std::vector connectionRates_; std::vector B_avg_; bool changed_to_stopped_this_step_ = false; diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index 839826d51..149a37119 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -48,13 +48,14 @@ WellInterfaceFluidSystem:: WellInterfaceFluidSystem(const Well& well, const ParallelWellInfo& 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>& perf_data) - : WellInterfaceGeneric(well, parallel_well_info, time_step, + : WellInterfaceGeneric(well, parallel_well_info, time_step, param, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data) , rateConverter_(rate_converter) diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.hpp b/opm/simulators/wells/WellInterfaceFluidSystem.hpp index 1a63d7bb4..a3dba1191 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.hpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.hpp @@ -57,6 +57,7 @@ protected: public: using Scalar = typename FluidSystem::Scalar; + using ModelParameters = typename WellInterfaceGeneric::ModelParameters; int flowPhaseToModelPhaseIdx(const int phaseIdx) const; @@ -73,6 +74,7 @@ protected: WellInterfaceFluidSystem(const Well& well, const ParallelWellInfo& parallel_well_info, const int time_step, + const ModelParameters& param, const RateConverterType& rate_converter, const int pvtRegionIdx, const int num_components, diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index c71fb90dd..8503522e0 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -57,6 +57,7 @@ WellInterfaceGeneric:: WellInterfaceGeneric(const Well& well, const ParallelWellInfo& 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& 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(max_inj_mult.size(), 1.); + const auto nperf = max_inj_mult.size(); + // initializing the (report?) step specific multipliers and damping factors to be 1.0 + this->inj_multiplier_ = std::vector(nperf, 1.); + this->inj_multiplier_previter_ = std::vector(nperf, 1.); + this->inj_multiplier_damp_factor_ = std::vector(nperf, 1.); } template @@ -243,13 +248,16 @@ template Scalar WellInterfaceGeneric:: getInjMult(const int perf, const Scalar bhp, - const Scalar perf_pres) const + const Scalar perf_pres, + DeferredLogger& dlogger) const { assert(!this->isProducer()); Scalar 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()) || @@ -265,15 +273,40 @@ 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) || + (multiplier < prev_multiplier && prev_multiplier > prev2_multiplier); + + 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) { + multiplier = std::max(multiplier, this->prev_inj_multiplier_[perf_ecl_index]); + } + + this->inj_multiplier_[perf_ecl_index] = multiplier; + this->inj_multiplier_previter_[perf_ecl_index] = 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) { - multiplier = std::max(multiplier, this->prev_inj_multiplier_[perf_ecl_index]); - } - - this->inj_multiplier_[perf_ecl_index] = multiplier; return multiplier; } diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index ce484b96b..e0956d70b 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -25,6 +25,7 @@ #define OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED #include +#include #include #include @@ -50,9 +51,12 @@ class Schedule; template class WellInterfaceGeneric { public: + using ModelParameters = BlackoilModelParameters; + WellInterfaceGeneric(const Well& well, const ParallelWellInfo& 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; @@ -274,6 +278,7 @@ protected: const ParallelWellInfo& 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 +360,11 @@ protected: // which intends to keep the fracturing open std::vector prev_inj_multiplier_; + // WINJMULT multipliers for previous iteration (used for oscillation detection) + mutable std::vector inj_multiplier_previter_; + // WINJMULT dampening factors (used in case of oscillations) + mutable std::vector inj_multiplier_damp_factor_; + // the multiplier due to injection filtration cake std::vector inj_fc_multiplier_; diff --git a/opm/simulators/wells/WellInterfaceIndices.cpp b/opm/simulators/wells/WellInterfaceIndices.cpp index a07dbdc1a..1150bdc18 100644 --- a/opm/simulators/wells/WellInterfaceIndices.cpp +++ b/opm/simulators/wells/WellInterfaceIndices.cpp @@ -38,6 +38,7 @@ WellInterfaceIndices:: WellInterfaceIndices(const Well& well, const ParallelWellInfo& parallel_well_info, const int time_step, + const ModelParameters& param, const typename WellInterfaceFluidSystem::RateConverterType& rate_converter, const int pvtRegionIdx, const int num_components, @@ -47,6 +48,7 @@ WellInterfaceIndices(const Well& well, : WellInterfaceFluidSystem(well, parallel_well_info, time_step, + param, rate_converter, pvtRegionIdx, num_components, diff --git a/opm/simulators/wells/WellInterfaceIndices.hpp b/opm/simulators/wells/WellInterfaceIndices.hpp index b339625c5..6ac6138e1 100644 --- a/opm/simulators/wells/WellInterfaceIndices.hpp +++ b/opm/simulators/wells/WellInterfaceIndices.hpp @@ -38,6 +38,7 @@ public: using WellInterfaceFluidSystem::Water; using Scalar = typename FluidSystem::Scalar; using Eval = DenseAd::Evaluation; + using ModelParameters = typename WellInterfaceFluidSystem::ModelParameters; int flowPhaseToModelCompIdx(const int phaseIdx) const; int modelCompIdxToFlowCompIdx(const unsigned compIdx) const; @@ -58,6 +59,7 @@ protected: WellInterfaceIndices(const Well& well, const ParallelWellInfo& parallel_well_info, const int time_step, + const ModelParameters& param, const typename WellInterfaceFluidSystem::RateConverterType& rate_converter, const int pvtRegionIdx, const int num_components, diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 92cac2d0a..b90256113 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -65,13 +65,13 @@ namespace Opm : WellInterfaceIndices(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; @@ -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; }