diff --git a/opm/simulators/flow/BlackoilModelParametersEbos.hpp b/opm/simulators/flow/BlackoilModelParametersEbos.hpp index d6f985d8bac..b0f00b39101 100644 --- a/opm/simulators/flow/BlackoilModelParametersEbos.hpp +++ b/opm/simulators/flow/BlackoilModelParametersEbos.hpp @@ -156,6 +156,7 @@ template struct MaxInnerIterWells { using type = UndefinedProperty; }; + template struct AlternativeWellRateInit { using type = UndefinedProperty; @@ -165,6 +166,24 @@ struct MaximumNumberOfWellSwitches { using type = UndefinedProperty; }; + +template +struct WellBhpScaling { + using type = double; + static constexpr type value = 1.0; +}; +template +struct WellRateScaling { + using type = double; + static constexpr type value = 1.0; +}; + +template +struct BhpControlScaling { + using type = double; + static constexpr type value = 1.0; +}; + template struct DbhpMaxRel { using type = GetPropType; @@ -495,6 +514,9 @@ namespace Opm EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking"); EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations"); EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control"); + EWOMS_REGISTER_PARAM(TypeTag, double, WellBhpScaling, "Scaling of the Bhp variable"); + EWOMS_REGISTER_PARAM(TypeTag, double, WellRateScaling, "Scaling of the Rate variable"); + EWOMS_REGISTER_PARAM(TypeTag, double, BhpControlScaling, "Scaling of bhp control equation"); } }; } // namespace Opm diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index d418a5db891..edee5589cda 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -350,7 +350,9 @@ namespace Opm { double gravity_{}; std::vector depth_{}; bool alternative_well_rate_init_{}; - + double bhp_scaling_; + double rate_scaling_; + double bhp_control_scaling_; std::unique_ptr rateConverter_{}; std::unique_ptr regionalAveragePressureCalculator_{}; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index b1c46ac9c98..9d50fc6df76 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -84,7 +84,11 @@ namespace Opm { extractLegacyDepth_(); gravity_ = ebosSimulator_.problem().gravity()[2]; - + // this parameters is registered in BlackoilModelParamtersEbos.hpp + // but no need of putting them in to the parameters. + bhp_scaling_ = EWOMS_GET_PARAM(TypeTag, double, WellBhpScaling); + rate_scaling_ = EWOMS_GET_PARAM(TypeTag, double, WellRateScaling); + bhp_control_scaling_ = EWOMS_GET_PARAM(TypeTag, double, BhpControlScaling); initial_step_ = true; // add the eWoms auxiliary module for the wells to the list @@ -106,6 +110,7 @@ namespace Opm { const bool well_opened_this_step = report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask); wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_, this->local_num_cells_, this->B_avg_, well_opened_this_step); + wellPtr->setScalings(bhp_scaling_,rate_scaling_,bhp_control_scaling_); } } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 669ce16a742..c9776cb698a 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -166,7 +166,9 @@ namespace Opm const SummaryState& summary_state, DeferredLogger& deferred_logger, double alq_value) const override; - + + void setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling); + double bhpControlScaling() const{return this->bhp_control_scaling_;} protected: int number_segments_; diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index c674450936f..3b242acfa85 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -214,7 +214,7 @@ checkConvergenceControlEq(const WellState& well_state, } } - const double well_control_residual = std::abs(resWell_[0][SPres]); + const double well_control_residual = std::abs(resWell_[0][SPres])/this->bhp_control_scaling_; const int dummy_component = -1; if (std::isnan(well_control_residual)) { report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()}); @@ -411,8 +411,8 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells, // update the segment pressure { const int sign = dwells[seg][SPres] > 0.? 1 : -1; - const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); - primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5); + const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change/this->bhp_scaling_); + primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5/this->bhp_scaling_); } // update the total rate // TODO: should we have a limitation of the total rate change? @@ -456,7 +456,7 @@ updatePrimaryVariables(const WellState& well_state) const // calculate the total rate for each segment double total_seg_rate = 0.0; // the segment pressure - primary_variables_[seg][SPres] = segment_pressure[seg]; + primary_variables_[seg][SPres] = segment_pressure[seg]/this->bhp_scaling_; // TODO: under what kind of circustances, the following will be wrong? // the definition of g makes the gas phase is always the last phase for (int p = 0; p < baseif_.numPhases(); p++) { @@ -470,7 +470,7 @@ updatePrimaryVariables(const WellState& well_state) const total_seg_rate = std::min(total_seg_rate, 0.); } } - primary_variables_[seg][WQTotal] = total_seg_rate; + primary_variables_[seg][WQTotal] = total_seg_rate/this->rate_scaling_; if (std::abs(total_seg_rate) > 0.) { if (has_wfrac_variable) { const int water_pos = pu.phase_pos[Water]; @@ -607,18 +607,18 @@ getSegmentRateUpwinding(const int seg, if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx && phase == InjectorType::WATER) - return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + return primary_variables_evaluation_[seg][WQTotal]*this->rate_scaling_ / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx && phase == InjectorType::OIL) - return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + return primary_variables_evaluation_[seg][WQTotal]*this->rate_scaling_ / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx && phase == InjectorType::GAS) - return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); + return primary_variables_evaluation_[seg][WQTotal]*this->rate_scaling_ / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx)); return 0.0; } @@ -627,7 +627,7 @@ getSegmentRateUpwinding(const int seg, assert(segment_rate.derivative(SPres + Indices::numEq) == 0.); - return segment_rate; + return segment_rate*this->rate_scaling_; } template @@ -821,7 +821,7 @@ typename MultisegmentWellEval::EvalWell MultisegmentWellEval:: getSegmentPressure(const int seg) const { - return primary_variables_evaluation_[seg][SPres]; + return primary_variables_evaluation_[seg][SPres]*this->bhp_scaling_; } template @@ -838,7 +838,7 @@ MultisegmentWellEval:: getSegmentRate(const int seg, const int comp_idx) const { - return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx); + return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx)*this->rate_scaling_; } template @@ -1543,9 +1543,9 @@ updateWellStateFromPrimaryVariables(WellState& well_state, } // update the segment pressure - segment_pressure[seg] = primary_variables_[seg][SPres]; + segment_pressure[seg] = primary_variables_[seg][SPres]*this->bhp_scaling_; if (seg == 0) { // top segment - ws.bhp = segment_pressure[seg]; + ws.bhp = segment_pressure[seg]*this->bhp_scaling_; } } updateThp(well_state, rho, deferred_logger); diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp index c7333e6e779..9cc9bbf708c 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.cpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -47,6 +47,9 @@ MultisegmentWellGeneric(WellInterfaceGeneric& baseif) , segment_inlets_(numberOfSegments()) , segment_depth_diffs_(numberOfSegments(), 0.0) , perforation_segment_depth_diffs_(baseif_.numPerfs(), 0.0) + , bhp_scaling_(1.0) + , rate_scaling_(1.0) + , bhp_control_scaling_(1.0) { // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. // for other facilities needed but not available from parser, we need to process them here diff --git a/opm/simulators/wells/MultisegmentWellGeneric.hpp b/opm/simulators/wells/MultisegmentWellGeneric.hpp index 9c8d2b47e32..4bb9f8c811c 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.hpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.hpp @@ -126,6 +126,9 @@ class MultisegmentWellGeneric // or in another way, the depth difference between the perforation and // the segment the perforation belongs to std::vector perforation_segment_depth_diffs_; + double bhp_scaling_; + double rate_scaling_; + double bhp_control_scaling_; }; } diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 0643073e6b0..1bd95f60aea 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -2028,6 +2028,12 @@ namespace Opm } + template + void MultisegmentWell::setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling){ + this->bhp_scaling_ = bhp_scaling; + this->rate_scaling_ = rate_scaling; + this->bhp_control_scaling_ = bhp_control_scaling; + } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 90d23a57d8a..1a174876166 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -249,7 +249,8 @@ namespace Opm double* connII, DeferredLogger& deferred_logger) const; - + void setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling); + double bhpControlScaling() const {return this->bhp_control_scaling_;} protected: bool regularize_; diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 25961913df8..318f4af6f1c 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -216,9 +216,9 @@ getQs(const int comp_idx) const // "Multi phase injectors are not supported, requested for well " + name()); break; } - return inj_frac * primary_variables_evaluation_[WQTotal]; + return inj_frac * this->getWQTotal();//primary_variables_evaluation_[WQTotal]*this->rate_scaling_; } else { // producers - return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx); + return this->getWQTotal() * wellVolumeFractionScaled(comp_idx); } } @@ -277,13 +277,13 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log if (baseif_.isInjector()) { switch (baseif_.wellEcl().injectorType()) { case InjectorType::WATER: - primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]]; + primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]]/this->rate_scaling_; break; case InjectorType::GAS: - primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]]; + primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]]/this->rate_scaling_; break; case InjectorType::OIL: - primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]]; + primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]]/this->rate_scaling_; break; case InjectorType::MULTI: // Not supported. @@ -292,7 +292,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log break; } } else { - primary_variables_[WQTotal] = total_well_rate; + primary_variables_[WQTotal] = total_well_rate/this->rate_scaling_; } if (std::abs(total_well_rate) > 0.) { @@ -353,7 +353,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log // BHP - primary_variables_[Bhp] = ws.bhp; + primary_variables_[Bhp] = ws.bhp/this->bhp_scaling_; } template @@ -665,12 +665,12 @@ updateWellStateFromPrimaryVariables(WellState& well_state, } auto& ws = well_state.well(baseif_.indexOfWell()); - ws.bhp = primary_variables_[Bhp]; + ws.bhp = primary_variables_[Bhp]*this->bhp_scaling_; // calculate the phase rates based on the primary variables // for producers, this is not a problem, while not sure for injectors here if (baseif_.isProducer()) { - const double g_total = primary_variables_[WQTotal]; + const double g_total = primary_variables_[WQTotal]*this->rate_scaling_; for (int p = 0; p < baseif_.numPhases(); ++p) { ws.surface_rates[p] = g_total * F[p]; } @@ -680,13 +680,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state, } switch (baseif_.wellEcl().injectorType()) { case InjectorType::WATER: - ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal]; + ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal]*this->rate_scaling_; break; case InjectorType::GAS: - ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal]; + ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal]*this->rate_scaling_; break; case InjectorType::OIL: - ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal]; + ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal]*this->rate_scaling_; break; case InjectorType::MULTI: // Not supported. @@ -771,9 +771,9 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells, // updating the bottom hole pressure { const int sign1 = dwells[0][Bhp] > 0 ? 1: -1; - const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit); + const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit);//NB dBHPLims ser param_.dbhp_max_rel_ // 1e5 to make sure bhp will not be below 1bar - primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5); + primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5/this->bhp_scaling_); } } diff --git a/opm/simulators/wells/StandardWellEval.hpp b/opm/simulators/wells/StandardWellEval.hpp index 224835abe75..5b5dc88271e 100644 --- a/opm/simulators/wells/StandardWellEval.hpp +++ b/opm/simulators/wells/StandardWellEval.hpp @@ -108,14 +108,14 @@ class StandardWellEval : public StandardWellGeneric void initPrimaryVariablesEvaluation() const; - const EvalWell& getBhp() const + const EvalWell getBhp() const { - return primary_variables_evaluation_[Bhp]; + return primary_variables_evaluation_[Bhp]*this->bhp_scaling_; } - const EvalWell& getWQTotal() const + const EvalWell getWQTotal() const { - return primary_variables_evaluation_[WQTotal]; + return primary_variables_evaluation_[WQTotal]*this->rate_scaling_; } EvalWell extendEval(const Eval& in) const; diff --git a/opm/simulators/wells/StandardWellGeneric.cpp b/opm/simulators/wells/StandardWellGeneric.cpp index a6de18d3241..95adb4891af 100644 --- a/opm/simulators/wells/StandardWellGeneric.cpp +++ b/opm/simulators/wells/StandardWellGeneric.cpp @@ -54,6 +54,10 @@ StandardWellGeneric(int Bhp, , perf_pressure_diffs_(baseif_.numPerfs()) , parallelB_(duneB_, baseif_.parallelWellInfo()) , Bhp_(Bhp) + , bhp_scaling_(1.0) + , rate_scaling_(1.0) + , bhp_control_scaling_(1.0) + { duneB_.setBuildMode(OffDiagMatWell::row_wise); duneC_.setBuildMode(OffDiagMatWell::row_wise); @@ -68,12 +72,13 @@ relaxationFactorRate(const std::vector& primary_variables, const BVectorWell& dwells) { double relaxation_factor = 1.0; - static constexpr int WQTotal = 0; + static constexpr int WQTotal = 0;//NB this may be inconsiten with other definitions // For injector, we only check the total rates to avoid sign change of rates - const double original_total_rate = primary_variables[WQTotal]; + // we only need scaled rates + const double original_total_rate = primary_variables[WQTotal];//*this->rate_scaling_; const double newton_update = dwells[0][WQTotal]; - const double possible_update_total_rate = primary_variables[WQTotal] - newton_update; + const double possible_update_total_rate = (primary_variables[WQTotal] - newton_update);;//*this->rate_scaling_; // 0.8 here is a experimental value, which remains to be optimized // if the original rate is zero or possible_update_total_rate is zero, relaxation_factor will @@ -466,7 +471,7 @@ checkConvergenceControlEq(const WellState& well_state, } } - const double well_control_residual = std::abs(this->resWell_[0][Bhp_]); + const double well_control_residual = std::abs(this->resWell_[0][Bhp_])/this->bhp_control_scaling_; const int dummy_component = -1; if (std::isnan(well_control_residual)) { report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()}); diff --git a/opm/simulators/wells/StandardWellGeneric.hpp b/opm/simulators/wells/StandardWellGeneric.hpp index d943d4234d1..6d10cfac586 100644 --- a/opm/simulators/wells/StandardWellGeneric.hpp +++ b/opm/simulators/wells/StandardWellGeneric.hpp @@ -32,7 +32,6 @@ #include #include - namespace Opm { @@ -137,6 +136,11 @@ class StandardWellGeneric private: int Bhp_; // index of Bhp +protected: + double bhp_scaling_; + double rate_scaling_; + double bhp_control_scaling_; + }; } diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index b1f9f6fad21..1ccc597f6c5 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -52,7 +52,7 @@ namespace Opm - + template void @@ -2081,7 +2081,12 @@ namespace Opm } } - + template + void StandardWell::setScalings(double bhp_scaling, double rate_scaling,double bhp_control_scaling){ + this->bhp_scaling_ = bhp_scaling; + this->rate_scaling_ = rate_scaling; + this->bhp_control_scaling_ = bhp_control_scaling; + } template diff --git a/opm/simulators/wells/WellInterfaceEval.cpp b/opm/simulators/wells/WellInterfaceEval.cpp index 0117fc8a85d..255e226ca28 100644 --- a/opm/simulators/wells/WellInterfaceEval.cpp +++ b/opm/simulators/wells/WellInterfaceEval.cpp @@ -47,7 +47,7 @@ namespace Opm template WellInterfaceEval:: WellInterfaceEval(const WellInterfaceFluidSystem& baseif) - : baseif_(baseif) + : baseif_(baseif) {} template @@ -342,11 +342,11 @@ assembleControlEqProd_(const WellState& well_state, break; } case Well::ProducerCMode::BHP: { - control_eq = bhp - controls.bhp_limit; + control_eq = (bhp - controls.bhp_limit)*this->bhpControlScaling(); break; } case Well::ProducerCMode::THP: { - control_eq = bhp - bhp_from_thp(); + control_eq = (bhp - bhp_from_thp())*this->bhpControlScaling(); break; } case Well::ProducerCMode::GRUP: { diff --git a/opm/simulators/wells/WellInterfaceEval.hpp b/opm/simulators/wells/WellInterfaceEval.hpp index 754e6567bd7..cffdfc01bc2 100644 --- a/opm/simulators/wells/WellInterfaceEval.hpp +++ b/opm/simulators/wells/WellInterfaceEval.hpp @@ -154,7 +154,7 @@ class WellInterfaceEval { const std::function& bhp_from_thp, EvalWell& control_eq, DeferredLogger& deferred_logger) const; - + virtual double bhpControlScaling() const = 0; protected: WellInterfaceEval(const WellInterfaceFluidSystem& baseif); diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 368693ad631..3c855e2949b 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -62,6 +62,8 @@ class WellInterfaceGeneric { const int index_of_well, const std::vector& perf_data); + virtual void setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling) = 0; + /// \brief Get the perforations of the well const std::vector& perforationData() const; @@ -360,6 +362,7 @@ class WellInterfaceGeneric { std::vector< std::string> well_control_log_; bool changed_to_open_this_step_ = true; + }; }