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); } }