Merge pull request #1045 from joakim-hove/guide-rate-nnn

Guide rate nnn

opm/parser/eclipse/EclipseState/Schedule/Group/Group2.hpp

@@ -89,13 +89,14 @@ enum class GuideRateTarget {
     WAT = 1,
     GAS = 2,
     LIQ = 3,
-    COMB = 4,
-    WGA =  5,
-    CVAL = 6,
-    INJV = 7,
-    POTN = 8,
-    FORM = 9,
-    NO_GUIDE_RATE = 10
+    RES = 4,
+    COMB = 5,
+    WGA =  6,
+    CVAL = 7,
+    INJV = 8,
+    POTN = 9,
+    FORM = 10,
+    NO_GUIDE_RATE = 11
 };
 static GuideRateTarget GuideRateTargetFromString( const std::string& stringValue );
 

opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp

@@ -25,16 +25,21 @@
 #include <ctime>
 #include <unordered_map>
 
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well2.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group/Group2.hpp>
+
 namespace Opm {
 
 class Schedule;
 class GuideRate {
+private:
 
 struct GuideRateValue {
     GuideRateValue() = default;
-    GuideRateValue(double t, double v):
+    GuideRateValue(double t, double v, GuideRateModel::Target tg):
         sim_time(t),
-        value(v)
+        value(v),
+        target(tg)
     {}
 
     bool operator==(const GuideRateValue& other) const {
@@ -46,23 +51,45 @@ struct GuideRateValue {
         return !(*this == other);
     }
 
+    double eval(GuideRateModel::Target target_arg) const;
+
     double sim_time;
     double value;
+    GuideRateModel::Target target;
 };
 
+struct Potential {
+    Potential() = default;
+    Potential(double op, double gp, double wp) :
+        oil_pot(op),
+        gas_pot(gp),
+        wat_pot(wp)
+    {}
+
+
+    double eval(Group2::GuideRateTarget target) const;
+    double eval(Well2::GuideRateTarget target) const;
+
+
+    double oil_pot;
+    double gas_pot;
+    double wat_pot;
+};
 
 public:
     GuideRate(const Schedule& schedule);
-    double update(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot);
+    void   compute(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot);
+    double get(const std::string& well, Well2::GuideRateTarget target) const;
+    double get(const std::string& group, Group2::GuideRateTarget target) const;
 private:
-    void well_update(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot);
-    void group_update(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot);
-    double get(const std::string& wgname) const;
+    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_group_pot() const;
     double eval_group_resvinj() const;
 
     std::unordered_map<std::string,GuideRateValue> values;
+    std::unordered_map<std::string,Potential> potentials;
     const Schedule& schedule;
 };
 

opm/parser/eclipse/EclipseState/Schedule/Group/GuideRateConfig.hpp

@@ -22,6 +22,7 @@
 
 #include <string>
 #include <unordered_map>
+#include <memory>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRateModel.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Group/Group2.hpp>
@@ -45,6 +46,7 @@ struct Group {
 };
 
     const GuideRateModel& model() const;
+    bool has_model() const;
     bool update_model(const GuideRateModel& model);
     void update_well(const Well2& well);
     void update_group(const Group2& group);
@@ -53,7 +55,7 @@ struct Group {
     bool has_well(const std::string& well) const;
     bool has_group(const std::string& group) const;
 private:
-    GuideRateModel m_model;
+    std::shared_ptr<GuideRateModel> m_model;
     std::unordered_map<std::string, Well> wells;
     std::unordered_map<std::string, Group> groups;
 };

opm/parser/eclipse/EclipseState/Schedule/Group/GuideRateModel.hpp

@@ -21,6 +21,8 @@
 #define GUIDE_RATE_MODEL_HPP
 
 #include <opm/parser/eclipse/Deck/UDAValue.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group/Group2.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well2.hpp>
 
 namespace Opm {
 
@@ -31,9 +33,10 @@ public:
         OIL = 0,
         LIQ = 1,
         GAS = 2,
-        RES = 3,
-        COMB = 4,
-        NONE = 5
+        WAT = 3,
+        RES = 4,
+        COMB = 5,
+        NONE = 6
     };
 
     static Target TargetFromString(const std::string& s);
@@ -59,7 +62,11 @@ public:
     bool operator==(const GuideRateModel& other) const;
     bool operator!=(const GuideRateModel& other) const;
     Target target() const;
+
+    static Target convert_target(Group2::GuideRateTarget group_target);
+    static Target convert_target(Well2::GuideRateTarget well_target);
 private:
+    double pot(double oil_pot, double gas_pot, double wat_pot) const;
     /*
       Unfortunately the default values will give a GuideRateModel which can not
       be evaluated, due to a division by zero problem.

opm/parser/eclipse/EclipseState/Schedule/Well/Well2.hpp

@@ -157,7 +157,7 @@ public:
         double scale_factor;
     };
 
-    
+
     struct InjectionControls {
     public:
         InjectionControls(int controls_arg) :

src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.cpp

@@ -22,74 +22,158 @@
 
 namespace Opm {
 
+double GuideRate::GuideRateValue::eval(GuideRateModel::Target target_arg) const
+{
+    if (target_arg == this->target)
+        return this->value;
+    else
+        throw std::logic_error("Don't know how to convert .... ");
+}
+
+double GuideRate::Potential::eval(Well2::GuideRateTarget target) const {
+    if (target == Well2::GuideRateTarget::OIL)
+        return this->oil_pot;
+
+    if (target == Well2::GuideRateTarget::GAS)
+        return this->gas_pot;
+
+    if (target == Well2::GuideRateTarget::LIQ)
+        return this->oil_pot + this->wat_pot;
+
+    if (target == Well2::GuideRateTarget::WAT)
+        return this->wat_pot;
+
+    throw std::logic_error("Don't know how to convert .... ");
+}
+
+double GuideRate::Potential::eval(Group2::GuideRateTarget target) const {
+    if (target == Group2::GuideRateTarget::OIL)
+        return this->oil_pot;
+
+    if (target == Group2::GuideRateTarget::GAS)
+        return this->gas_pot;
+
+    if (target == Group2::GuideRateTarget::LIQ)
+        return this->oil_pot + this->wat_pot;
+
+    if (target == Group2::GuideRateTarget::WAT)
+        return this->wat_pot;
+
+    throw std::logic_error("Don't know how to convert .... ");
+}
+
+
 GuideRate::GuideRate(const Schedule& schedule_arg) :
     schedule(schedule_arg)
 {}
 
 
-double GuideRate::get(const std::string& wgname) const {
-    const auto& value = this->values.at(wgname);
-    return value.value;
+
+double GuideRate::get(const std::string& well, Well2::GuideRateTarget target) const {
+    const auto iter = this->values.find(well);
+    if (iter != this->values.end()) {
+        const auto& value = iter->second;
+        auto model_target = GuideRateModel::convert_target(target);
+        return value.eval(model_target);
+    } else {
+        const auto& pot = this->potentials.at(well);
+        return pot.eval(target);
+    }
+}
+
+double GuideRate::get(const std::string& group, Group2::GuideRateTarget target) const {
+    const auto iter = this->values.find(group);
+     if (iter != this->values.end()) {
+        auto model_target = GuideRateModel::convert_target(target);
+        const auto& value = this->values.at(group);
+        return value.eval(model_target);
+    } else {
+        const auto& pot = this->potentials.at(group);
+        return pot.eval(target);
+    }
 }
 
 
-double GuideRate::update(const std::string& wgname, size_t report_step, double sim_time, double oil_pot, double gas_pot, double wat_pot) {
+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);
+    this->potentials[wgname] = Potential{oil_pot, gas_pot, wat_pot};
 
-    if (config.has_well(wgname))
-        this->well_update(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot);
-    else if (config.has_group(wgname)) {
-        this->group_update(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot);
-    } else
-        throw std::out_of_range("No such well/group: ");
+    if (config.has_group(wgname))
+        this->group_compute(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot);
+    else
+        this->well_compute(wgname, report_step, sim_time, oil_pot, gas_pot, wat_pot);
 
-    return this->get(wgname);
 }
 
 
-void GuideRate::group_update(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);
-    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 (iter != this->values.end()) {
-        const auto& grv = iter->second;
+    if (group.guide_rate > 0) {
+        auto model_target = GuideRateModel::convert_target(group.target);
+        this->values[wgname] = GuideRateValue( 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 (iter != this->values.end()) {
+            const auto& grv = iter->second;
+            if (group.target == Group2::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");
+
+                auto time_diff = sim_time - grv.sim_time;
+                if (config.model().update_delay() > time_diff)
+                    return;
+            }
+        }
+
+
+        if (group.target == Group2::GuideRateTarget::INJV)
+            throw std::logic_error("Group guide rate mode: INJV not implemented");
+
+        if (group.target == Group2::GuideRateTarget::POTN)
+            throw std::logic_error("Group guide rate mode: POTN not implemented");
+
         if (group.target == Group2::GuideRateTarget::FORM) {
-            auto time_diff = sim_time - grv.sim_time;
-            if (config.model().update_delay() > time_diff)
-                return;
+            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 (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()};
         }
     }
-
-
-    double guide_rate = group.guide_rate;
-
-    if (group.guide_rate == 0 || group.target == Group2::GuideRateTarget::POTN)
-        guide_rate = this->eval_group_pot();
-
-    if (group.target == Group2::GuideRateTarget::INJV)
-        guide_rate = this->eval_group_resvinj();
-
-    if (group.target == Group2::GuideRateTarget::FORM) {
-        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));
-    }
-
-    this->values[wgname] = GuideRateValue{sim_time, guide_rate};
 }
 
 
-void GuideRate::well_update(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);
-    const auto& well = config.well(wgname);
 
-    if (well.guide_rate > 0)
-        this->values[wgname] = GuideRateValue( sim_time, well.guide_rate );
-    else {
+    // guide rates spesified with WGRUPCON
+    if (config.has_well(wgname)) {
+        const auto& well = config.well(wgname);
+        if (well.guide_rate > 0) {
+            auto model_target = GuideRateModel::convert_target(well.target);
+            this->values[wgname] = GuideRateValue( sim_time, well.guide_rate, model_target );
+        }
+    } else if (config.has_model()) { // GUIDERAT
+        // only look for wells not groups
+        if (!this->schedule.hasWell(wgname, report_step))
+            return;
+
+        const auto& well = this->schedule.getWell2(wgname, report_step);
+
+        // GUIDERAT does not apply to injectors
+        if (well.isInjector())
+            return;
+
         auto iter = this->values.find(wgname);
 
         if (iter != this->values.end()) {
@@ -105,8 +189,9 @@ void GuideRate::well_update(const std::string& wgname, size_t report_step, doubl
         else
             guide_rate = this->eval_form(config.model(),  oil_pot,  gas_pot,  wat_pot, std::addressof(iter->second));
 
-        this->values[wgname] = GuideRateValue{sim_time, guide_rate};
+        this->values[wgname] = GuideRateValue{sim_time, guide_rate, config.model().target()};
     }
+    // If neither WGRUPCON nor GUIDERAT is spesified potentials are used
 }
 
 double GuideRate::eval_form(const GuideRateModel& model, double oil_pot, double gas_pot, double wat_pot, const GuideRateValue * prev) const {

src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRateConfig.cpp

@@ -25,13 +25,23 @@ namespace Opm {
 
 
 const GuideRateModel& GuideRateConfig::model() const {
-    return this->m_model;
+    if (this->m_model)
+        return *this->m_model;
+    else
+        throw std::logic_error("Tried to dereference empty GuideRateModel");
 }
 
 
+bool GuideRateConfig::has_model() const {
+    if (this->m_model)
+        return true;
+    else
+        return false;
+}
+
 bool GuideRateConfig::update_model(const GuideRateModel& new_model) {
-    if (this->m_model != new_model) {
-        this->m_model = new_model;
+    if (!this->m_model || *(this->m_model) != new_model) {
+        this->m_model.reset( new GuideRateModel(new_model) );
         return true;
     }
     return false;

src/opm/parser/eclipse/EclipseState/Schedule/Group/GuideRateModel.cpp

@@ -68,32 +68,52 @@ GuideRateModel::GuideRateModel(double time_interval_arg,
         throw std::logic_error("Sorry - the 'COMB' mode is not supported");
 }
 
+double GuideRateModel::pot(double oil_pot, double gas_pot, double wat_pot) const {
+    switch (this->m_target) {
+    case Target::OIL:
+        return oil_pot;
 
+    case Target::LIQ:
+        return oil_pot + wat_pot;
+
+    case Target::GAS:
+        return gas_pot;
+
+    case Target::COMB:
+        throw std::logic_error("Not implemented - don't have a clue?");
+
+    case Target::RES:
+        throw std::logic_error("Not implemented - don't have a clue?");
+
+    default:
+        throw std::logic_error("Hmmm - should not be here?");
+    }
+}
 
 
 double GuideRateModel::eval(double oil_pot, double gas_pot, double wat_pot) const {
     if (this->default_model)
         throw std::invalid_argument("The default GuideRateModel can not be evaluated - must enter GUIDERAT information explicitly.");
 
-    double pot;
+    double pot = this->pot(oil_pot, gas_pot, wat_pot);
+    if (pot == 0)
+        return 0;
+
+
     double R1;
     double R2;
-
     switch (this->m_target) {
     case Target::OIL:
-        pot = oil_pot;
         R1 = wat_pot / oil_pot;
         R2 = gas_pot / oil_pot;
         break;
 
     case Target::LIQ:
-        pot = oil_pot + wat_pot;;
         R1 = oil_pot / (oil_pot + wat_pot);
         R2 = wat_pot / (oil_pot + wat_pot);
         break;
 
     case Target::GAS:
-        pot = gas_pot;
         R1 = wat_pot / gas_pot;
         R2 = oil_pot / gas_pot;
         break;
@@ -200,5 +220,29 @@ bool GuideRateModel::updateLINCOM(const UDAValue& , const UDAValue& , const UDAV
 }
 
 
+GuideRateModel::Target GuideRateModel::convert_target(Well2::GuideRateTarget well_target) {
+    if (well_target == Well2::GuideRateTarget::OIL)
+        return Target::OIL;
 
+    if (well_target == Well2::GuideRateTarget::GAS)
+        return Target::GAS;
+
+    if (well_target == Well2::GuideRateTarget::LIQ)
+        return Target::LIQ;
+
+    throw std::logic_error("Can not convert this .... ");
+}
+
+GuideRateModel::Target GuideRateModel::convert_target(Group2::GuideRateTarget group_target) {
+    if (group_target == Group2::GuideRateTarget::OIL)
+        return Target::OIL;
+
+    if (group_target == Group2::GuideRateTarget::GAS)
+        return Target::GAS;
+
+    if (group_target == Group2::GuideRateTarget::LIQ)
+        return Target::LIQ;
+
+    throw std::logic_error("Can not convert this .... ");
+}
 }

tests/parser/ScheduleTests.cpp

@@ -3565,7 +3565,8 @@ DATES             -- 4
     double oil_pot = 1;
     double gas_pot = 1;
     double wat_pot = 1;
-    BOOST_CHECK_THROW(gr.update("XYZ",1, 1.0, oil_pot, gas_pot, wat_pot), std::out_of_range);
+
+    gr.compute("XYZ",1, 1.0, oil_pot, gas_pot, wat_pot);
 }