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.
This commit is contained in:
Bård Skaflestad 2020-08-17 23:37:27 +02:00
parent 3080b9c7f6
commit 7649df53a2
2 changed files with 175 additions and 104 deletions

View File

@ -22,6 +22,8 @@
#include <cstddef> #include <cstddef>
#include <ctime> #include <ctime>
#include <limits>
#include <memory>
#include <string> #include <string>
#include <unordered_map> #include <unordered_map>
@ -77,11 +79,15 @@ struct GuideRateValue {
return !(*this == other); return !(*this == other);
} }
double sim_time; double sim_time { std::numeric_limits<double>::lowest() };
double value; double value { std::numeric_limits<double>::lowest() };
GuideRateModel::Target target; GuideRateModel::Target target { GuideRateModel::Target::NONE };
}; };
struct GRValState {
GuideRateValue curr{};
GuideRateValue prev{};
};
public: public:
GuideRate(const Schedule& schedule); GuideRate(const Schedule& schedule);
@ -94,11 +100,16 @@ public:
private: 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 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); 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_pot() const;
double eval_group_resvinj() const; double eval_group_resvinj() const;
std::unordered_map<std::string, GuideRateValue> 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<GRValState>;
std::unordered_map<std::string, GRValPtr> values;
std::unordered_map<std::string, RateVector > potentials; std::unordered_map<std::string, RateVector > potentials;
const Schedule& schedule; const Schedule& schedule;
}; };

View File

@ -20,11 +20,18 @@
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp> #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp> #include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <algorithm>
#include <memory>
#include <stdexcept>
#include <string>
#include <utility>
#include <stddef.h>
namespace Opm { namespace Opm {
double GuideRate::RateVector::eval(Well::GuideRateTarget target) const double GuideRate::RateVector::eval(Well::GuideRateTarget target) const
{ {
return this->eval(GuideRateModel::convert_target(target)); return this->eval(GuideRateModel::convert_target(target));
@ -59,40 +66,41 @@ GuideRate::GuideRate(const Schedule& schedule_arg) :
schedule(schedule_arg) schedule(schedule_arg)
{} {}
double GuideRate::get(const std::string& well, Well::GuideRateTarget target, const RateVector& rates) const
{
double GuideRate::get(const std::string& well, Well::GuideRateTarget target, const RateVector& rates) const { return this->get(well, GuideRateModel::convert_target(target), rates);
auto model_target = GuideRateModel::convert_target(target);
return get(well, model_target, rates);
} }
double GuideRate::get(const std::string& group, Group::GuideRateTarget target, const RateVector& rates) const { 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); return this->get(group, GuideRateModel::convert_target(target), rates);
} }
double GuideRate::get(const std::string& name, GuideRateModel::Target model_target, const RateVector& rates) const { 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()) { using namespace unit;
const auto& value = iter->second; using prefix::micro;
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) auto iter = this->values.find(name);
return value.value; if (iter == this->values.end()) {
if (value_target_rate < 1e-6) return this->potentials.at(name).eval(model_target);
return value.value;
// 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);
} }
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 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; return this->values.count(name) > 0;
} }
void GuideRate::compute(const std::string& wgname,
size_t 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) { double sim_time,
const auto& config = this->schedule.guideRateConfig(report_step); double oil_pot,
double gas_pot,
double wat_pot)
{
this->potentials[wgname] = RateVector{oil_pot, gas_pot, 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); 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); this->well_compute(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot);
}
} }
void GuideRate::group_compute(const std::string& wgname,
void GuideRate::group_compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot) { 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& config = this->schedule.guideRateConfig(report_step);
const auto& group = config.group(wgname); 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); 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); auto iter = this->values.find(wgname);
// If the FORM mode is used we check if the last computation is recent enough; // If the FORM mode is used we check if the last computation is
// then we just return. // recent enough; then we just return.
if (iter != this->values.end()) { if (iter != this->values.end()) {
const auto& grv = iter->second;
if (group.target == Group::GuideRateTarget::FORM) { if (group.target == Group::GuideRateTarget::FORM) {
if (!config.has_model()) if (! config.has_model()) {
throw std::logic_error("When specifying GUIDERATE target FORM you must enter a guiderate model with the GUIDERAT keyword"); 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; const auto& grv = iter->second->curr;
if (config.model().update_delay() > time_diff) const auto time_diff = sim_time - grv.sim_time;
if (config.model().update_delay() > time_diff) {
return; return;
}
} }
} }
if (group.target == Group::GuideRateTarget::INJV) {
if (group.target == Group::GuideRateTarget::INJV)
throw std::logic_error("Group guide rate mode: INJV not implemented"); 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"); throw std::logic_error("Group guide rate mode: POTN not implemented");
}
if (group.target == Group::GuideRateTarget::FORM) { if (group.target == Group::GuideRateTarget::FORM) {
double guide_rate; if (! config.has_model()) {
if (!config.has_model()) throw std::logic_error {
throw std::logic_error("When specifying GUIDERATE target FORM you must enter a guiderate model with the GUIDERAT keyword"); "When specifying GUIDERATE target FORM you must "
"enter a guiderate model with the GUIDERAT keyword"
};
}
if (iter != this->values.end()) const auto guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot);
guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, std::addressof(iter->second)); this->assign_grvalue(wgname, config.model(), { sim_time, guide_rate, config.model().target() });
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()};
} }
} }
} }
void GuideRate::well_compute(const std::string& wgname,
void GuideRate::well_compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot) { 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& config = this->schedule.guideRateConfig(report_step);
// guide rates spesified with WGRUPCON // guide rates spesified with WGRUPCON
if (config.has_well(wgname)) { if (config.has_well(wgname)) {
const auto& well = config.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); 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 // only look for wells not groups
if (!this->schedule.hasWell(wgname, report_step)) if (! this->schedule.hasWell(wgname, report_step)) {
return; return;
}
const auto& well = this->schedule.getWell(wgname, report_step); const auto& well = this->schedule.getWell(wgname, report_step);
// GUIDERAT does not apply to injectors // GUIDERAT does not apply to injectors
if (well.isInjector()) if (well.isInjector()) {
return; 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; auto iter = this->values.find(wgname);
if (iter == this->values.end()) if (iter != this->values.end()) {
guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, nullptr); const auto& grv = iter->second->curr;
else const auto time_diff = sim_time - grv.sim_time;
guide_rate = this->eval_form(config.model(), oil_pot, gas_pot, wat_pot, std::addressof(iter->second)); 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 GuideRate::eval_form(const GuideRateModel& model, double oil_pot, double gas_pot, double wat_pot) const
double new_guide_rate = model.eval(oil_pot, gas_pot, wat_pot); {
if (!prev) return model.eval(oil_pot, gas_pot, wat_pot);
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_group_pot() const { double GuideRate::eval_group_pot() const
return 0; {
return 0.0;
} }
double GuideRate::eval_group_resvinj() const { double GuideRate::eval_group_resvinj() const
return 0; {
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<GRValState>();
}
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);
} }
} }