From 7649df53a2268a5b7525fab633c1191d012a25b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 17 Aug 2020 23:37:27 +0200 Subject: [PATCH] GuideRate: Revise Damping Factor Implementation This commit centralises the way we incorporate damping factors (item 10 of the GUIDERAT keyword) into the calculation of group/well guide rates. In particular, we create a structure that manages both the current and the previous (damped) guiderate values and ensures that the new guiderate value is GR_p = f*GR_p' + (1 - f)*GR_p^{n-1} with GR_p' denoting the "raw" guiderate value calculated directly from potential rates at the current timelevel (n). GR_p^{n-1} is the damped-and previously used-guiderate value from timelevel n-1. Finally 'f' denotes the damping factor. This is the same approach used previously, but with some small changes to exclude zero-valued guiderates. We furthermore remove one of the early returns in GuideRate::get(). There is no need to return the nominated phase's guiderate value if the model phase rate is very low and doing so produces incorrect water guiderates for the OPL5 well in the MOD4_GRP test case. --- .../EclipseState/Schedule/Group/GuideRate.hpp | 21 +- .../EclipseState/Schedule/Group/GuideRate.cpp | 258 +++++++++++------- 2 files changed, 175 insertions(+), 104 deletions(-) diff --git a/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp b/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp index be9005ee0..73f160e6d 100644 --- a/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp +++ b/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp @@ -22,6 +22,8 @@ #include #include +#include +#include #include #include @@ -77,11 +79,15 @@ struct GuideRateValue { return !(*this == other); } - double sim_time; - double value; - GuideRateModel::Target target; + double sim_time { std::numeric_limits::lowest() }; + double value { std::numeric_limits::lowest() }; + GuideRateModel::Target target { GuideRateModel::Target::NONE }; }; +struct GRValState { + GuideRateValue curr{}; + GuideRateValue prev{}; +}; public: GuideRate(const Schedule& schedule); @@ -94,11 +100,16 @@ public: private: void well_compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot); void group_compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot); - double eval_form(const GuideRateModel& model, double oil_pot, double gas_pot, double wat_pot, const GuideRateValue * prev) const; + double eval_form(const GuideRateModel& model, double oil_pot, double gas_pot, double wat_pot) const; double eval_group_pot() const; double eval_group_resvinj() const; - std::unordered_map values; + void assign_grvalue(const std::string& wgname, const GuideRateModel& model, GuideRateValue&& value); + double get_grvalue_result(const GRValState& gr) const; + + using GRValPtr = std::unique_ptr; + + std::unordered_map values; std::unordered_map potentials; const Schedule& schedule; }; diff --git a/src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.cpp b/src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.cpp index 7e31af0a0..8955cc2ef 100644 --- a/src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.cpp +++ b/src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.cpp @@ -20,11 +20,18 @@ #include #include +#include + +#include +#include +#include +#include +#include + +#include + namespace Opm { - - - double GuideRate::RateVector::eval(Well::GuideRateTarget target) const { return this->eval(GuideRateModel::convert_target(target)); @@ -59,40 +66,41 @@ GuideRate::GuideRate(const Schedule& schedule_arg) : schedule(schedule_arg) {} - - -double GuideRate::get(const std::string& well, Well::GuideRateTarget target, const RateVector& rates) const { - auto model_target = GuideRateModel::convert_target(target); - return get(well, model_target, rates); +double GuideRate::get(const std::string& well, Well::GuideRateTarget target, const RateVector& rates) const +{ + return this->get(well, GuideRateModel::convert_target(target), rates); } -double GuideRate::get(const std::string& group, Group::GuideRateTarget target, const RateVector& rates) const { - auto model_target = GuideRateModel::convert_target(target); - return get(group, model_target, rates); +double GuideRate::get(const std::string& group, Group::GuideRateTarget target, const RateVector& rates) const +{ + return this->get(group, GuideRateModel::convert_target(target), rates); } -double GuideRate::get(const std::string& name, GuideRateModel::Target model_target, const RateVector& rates) const { - const auto iter = this->values.find(name); - if (iter != this->values.end()) { - const auto& value = iter->second; - if (value.target == model_target) - return value.value; - else { - double model_target_rate = rates.eval(model_target); - double value_target_rate = rates.eval(value.target); - - if (model_target_rate < 1e-6) - return value.value; - if (value_target_rate < 1e-6) - return value.value; +double GuideRate::get(const std::string& name, GuideRateModel::Target model_target, const RateVector& rates) const +{ + using namespace unit; + using prefix::micro; - // scale with the current production ratio when the control target differs from the guide rate target. - return value.value * model_target_rate / value_target_rate; - } - } else { - const auto& pot = this->potentials.at(name); - return pot.eval(model_target); + auto iter = this->values.find(name); + if (iter == this->values.end()) { + return this->potentials.at(name).eval(model_target); } + + const auto& value = *iter->second; + const auto grvalue = this->get_grvalue_result(value); + if (value.curr.target == model_target) { + return grvalue; + } + + const auto value_target_rate = rates.eval(value.curr.target); + if (value_target_rate < 1.0*micro*liter/day) { + return grvalue; + } + + // Scale with the current production ratio when the control target + // differs from the guide rate target. + const auto scale = rates.eval(model_target) / value_target_rate; + return grvalue * scale; } bool GuideRate::has(const std::string& name) const @@ -100,127 +108,179 @@ bool GuideRate::has(const std::string& name) const return this->values.count(name) > 0; } - - -void GuideRate::compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot) { - const auto& config = this->schedule.guideRateConfig(report_step); +void GuideRate::compute(const std::string& wgname, + size_t report_step, + double sim_time, + double oil_pot, + double gas_pot, + double wat_pot) +{ this->potentials[wgname] = RateVector{oil_pot, gas_pot, wat_pot}; - if (config.has_group(wgname)) + const auto& config = this->schedule.guideRateConfig(report_step); + if (config.has_group(wgname)) { this->group_compute(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot); - else + } + else { this->well_compute(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot); - + } } - -void GuideRate::group_compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot) { +void GuideRate::group_compute(const std::string& wgname, + size_t report_step, + double sim_time, + double oil_pot, + double gas_pot, + double wat_pot) +{ const auto& config = this->schedule.guideRateConfig(report_step); const auto& group = config.group(wgname); - if (group.guide_rate > 0) { + if (group.guide_rate > 0.0) { auto model_target = GuideRateModel::convert_target(group.target); - this->values[wgname] = GuideRateValue( sim_time, group.guide_rate, model_target ); - } else { + + const auto& model = config.has_model() ? config.model() : GuideRateModel{}; + this->assign_grvalue(wgname, model, { sim_time, group.guide_rate, model_target }); + } + else { auto iter = this->values.find(wgname); - // If the FORM mode is used we check if the last computation is recent enough; - // then we just return. + // If the FORM mode is used we check if the last computation is + // recent enough; then we just return. if (iter != this->values.end()) { - const auto& grv = iter->second; if (group.target == Group::GuideRateTarget::FORM) { - if (!config.has_model()) - throw std::logic_error("When specifying GUIDERATE target FORM you must enter a guiderate model with the GUIDERAT keyword"); + if (! config.has_model()) { + throw std::logic_error { + "When specifying GUIDERATE target FORM you must " + "enter a guiderate model with the GUIDERAT keyword" + }; + } - auto time_diff = sim_time - grv.sim_time; - if (config.model().update_delay() > time_diff) + const auto& grv = iter->second->curr; + const auto time_diff = sim_time - grv.sim_time; + if (config.model().update_delay() > time_diff) { return; + } } } - - if (group.target == Group::GuideRateTarget::INJV) + if (group.target == Group::GuideRateTarget::INJV) { throw std::logic_error("Group guide rate mode: INJV not implemented"); + } - if (group.target == Group::GuideRateTarget::POTN) + if (group.target == Group::GuideRateTarget::POTN) { throw std::logic_error("Group guide rate mode: POTN not implemented"); + } if (group.target == Group::GuideRateTarget::FORM) { - double guide_rate; - if (!config.has_model()) - throw std::logic_error("When specifying GUIDERATE target FORM you must enter a guiderate model with the GUIDERAT keyword"); + if (! config.has_model()) { + throw std::logic_error { + "When specifying GUIDERATE target FORM you must " + "enter a guiderate model with the GUIDERAT keyword" + }; + } - if (iter != this->values.end()) - guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, std::addressof(iter->second)); - else - guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, nullptr); - - this->values[wgname] = GuideRateValue{sim_time, guide_rate, config.model().target()}; + const auto guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot); + this->assign_grvalue(wgname, config.model(), { sim_time, guide_rate, config.model().target() }); } } } - -void GuideRate::well_compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot) { +void GuideRate::well_compute(const std::string& wgname, + size_t report_step, + double sim_time, + double oil_pot, + double gas_pot, + double wat_pot) +{ const auto& config = this->schedule.guideRateConfig(report_step); // guide rates spesified with WGRUPCON if (config.has_well(wgname)) { const auto& well = config.well(wgname); - if (well.guide_rate > 0) { + if (well.guide_rate > 0.0) { auto model_target = GuideRateModel::convert_target(well.target); - this->values[wgname] = GuideRateValue( sim_time, well.guide_rate, model_target ); + + const auto& model = config.has_model() ? config.model() : GuideRateModel{}; + this->assign_grvalue(wgname, model, { sim_time, well.guide_rate, model_target }); } - } else if (config.has_model()) { // GUIDERAT + } + else if (config.has_model()) { // GUIDERAT // only look for wells not groups - if (!this->schedule.hasWell(wgname, report_step)) + if (! this->schedule.hasWell(wgname, report_step)) { return; + } const auto& well = this->schedule.getWell(wgname, report_step); // GUIDERAT does not apply to injectors - if (well.isInjector()) + if (well.isInjector()) { return; - - auto iter = this->values.find(wgname); - - if (iter != this->values.end()) { - const auto& grv = iter->second; - auto time_diff = sim_time - grv.sim_time; - if (config.model().update_delay() > time_diff) - return; } - double guide_rate; - if (iter == this->values.end()) - guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, nullptr); - else - guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, std::addressof(iter->second)); + auto iter = this->values.find(wgname); + if (iter != this->values.end()) { + const auto& grv = iter->second->curr; + const auto time_diff = sim_time - grv.sim_time; + if (config.model().update_delay() > time_diff) { + return; + } + } - this->values[wgname] = GuideRateValue{sim_time, guide_rate, config.model().target()}; + const auto guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot); + this->assign_grvalue(wgname, config.model(), { sim_time, guide_rate, config.model().target() }); } - // If neither WGRUPCON nor GUIDERAT is spesified potentials are used + // If neither WGRUPCON nor GUIDERAT is specified potentials are used } -double GuideRate::eval_form(const GuideRateModel& model, double oil_pot, double gas_pot, double wat_pot, const GuideRateValue * prev) const { - double new_guide_rate = model.eval(oil_pot, gas_pot, wat_pot); - if (!prev) - return new_guide_rate; - - if (new_guide_rate > prev->value && !model.allow_increase()) - new_guide_rate = prev->value; - - auto damping_factor = model.damping_factor(); - - return damping_factor * new_guide_rate + (1 - damping_factor) * prev->value; +double GuideRate::eval_form(const GuideRateModel& model, double oil_pot, double gas_pot, double wat_pot) const +{ + return model.eval(oil_pot, gas_pot, wat_pot); } -double GuideRate::eval_group_pot() const { - return 0; +double GuideRate::eval_group_pot() const +{ + return 0.0; } -double GuideRate::eval_group_resvinj() const { - return 0; +double GuideRate::eval_group_resvinj() const +{ + return 0.0; +} + +void GuideRate::assign_grvalue(const std::string& wgname, + const GuideRateModel& model, + GuideRateValue&& value) +{ + using std::swap; + + auto& v = this->values[wgname]; + if (v == nullptr) { + v = std::make_unique(); + } + + swap(v->prev, v->curr); + + v->curr = std::move(value); + + if ((v->prev.sim_time < 0.0) || ! (v->prev.value > 0.0)) { + // No previous non-zero guiderate exists. No further actions. + return; + } + + // Incorporate damping &c. + const auto new_guide_rate = model.allow_increase() + ? v->curr.value : std::min(v->curr.value, v->prev.value); + + const auto damping_factor = model.damping_factor(); + v->curr.value = damping_factor*new_guide_rate + (1 - damping_factor)*v->prev.value; +} + +double GuideRate::get_grvalue_result(const GRValState& gr) const +{ + return (gr.curr.sim_time < 0.0) + ? 0.0 + : std::max(gr.curr.value, 0.0); } }