Merge pull request #4669 from steink/explicit_vfp_fallback

Include fallback to explicit vfp lookup for tiny rates - testing

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -235,7 +235,7 @@ namespace Opm {
                 // update VFP properties
                 vfp_properties_ = std::make_unique<VFPProperties>(sched_state.vfpinj(),
                                                                   sched_state.vfpprod(),
-                                                                  this->prevWellState());
+                                                                  this->wellState());
                 this->initializeWellProdIndCalculators();
                 if (sched_state.events().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX)) {
                     this->runWellPIScaling(timeStepIdx, local_deferredLogger);
@@ -340,6 +340,12 @@ namespace Opm {
             }
         }
 
+        for (auto& well : well_container_) {
+            if (well->isVFPActive(local_deferredLogger)){
+                well->setPrevSurfaceRates(this->wellState(), this->prevWellState());
+            }
+        }
+
         // calculate the well potentials
         try {
             updateWellPotentials(reportStepIdx,

opm/simulators/wells/SingleWellState.cpp

@@ -44,6 +44,7 @@ SingleWellState::SingleWellState(const std::string& name_,
     , productivity_index(pu_.num_phases)
     , surface_rates(pu_.num_phases)
     , reservoir_rates(pu_.num_phases)
+    , prev_surface_rates(pu_.num_phases)
     , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases)
     , trivial_target(false)
 {
@@ -84,6 +85,7 @@ void SingleWellState::shut() {
     this->thp = 0;
     this->status = Well::Status::SHUT;
     std::fill(this->surface_rates.begin(), this->surface_rates.end(), 0);
+    std::fill(this->prev_surface_rates.begin(), this->prev_surface_rates.end(), 0);
     std::fill(this->reservoir_rates.begin(), this->reservoir_rates.end(), 0);
     std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0);
 
@@ -296,6 +298,7 @@ bool SingleWellState::operator==(const SingleWellState& rhs) const
            this->productivity_index == rhs.productivity_index &&
            this->surface_rates == rhs.surface_rates &&
            this->reservoir_rates == rhs.reservoir_rates &&
+           this->prev_surface_rates == rhs.prev_surface_rates &&
            this->trivial_target == rhs.trivial_target &&
            this->segments == rhs.segments &&
            this->events == rhs.events &&

opm/simulators/wells/SingleWellState.hpp

@@ -63,6 +63,7 @@ public:
         serializer(productivity_index);
         serializer(surface_rates);
         serializer(reservoir_rates);
+        serializer(prev_surface_rates);
         serializer(trivial_target);
         serializer(segments);
         serializer(events);
@@ -95,6 +96,7 @@ public:
     std::vector<double> productivity_index;
     std::vector<double> surface_rates;
     std::vector<double> reservoir_rates;
+    std::vector<double> prev_surface_rates;
     PerfData perf_data;
     bool trivial_target;
     SegmentState segments;

opm/simulators/wells/VFPHelpers.cpp

@@ -348,7 +348,7 @@ VFPEvaluation bhp(const VFPProdTable& table,
     double flo = detail::getFlo(table, aqua, liquid, vapour);
     double wfr = detail::getWFR(table, aqua, liquid, vapour);
     double gfr = detail::getGFR(table, aqua, liquid, vapour);
-    if (use_vfpexplicit) {
+    if (use_vfpexplicit || -flo < table.getFloAxis().front()) {
         wfr = explicit_wfr;
         gfr = explicit_gfr;
     }

opm/simulators/wells/VFPProdProperties.cpp

@@ -161,7 +161,7 @@ EvalWell VFPProdProperties::bhp(const int table_id,
     EvalWell flo = detail::getFlo(table, aqua, liquid, vapour);
     EvalWell wfr = detail::getWFR(table, aqua, liquid, vapour);
     EvalWell gfr = detail::getGFR(table, aqua, liquid, vapour);
-    if (use_expvfp) {
+    if (use_expvfp || -flo.value() < table.getFloAxis().front()) {
         wfr = explicit_wfr;
         gfr = explicit_gfr;
     }

opm/simulators/wells/VFPProperties.hpp

@@ -72,7 +72,7 @@ public:
     }
 
     double getExplicitWFR(const int table_id, const size_t well_index) const {
-        const auto& rates = well_state_.well(well_index).surface_rates;
+        const auto& rates = well_state_.well(well_index).prev_surface_rates;
         const auto& pu = well_state_.phaseUsage();
         const auto& aqua = pu.phase_used[BlackoilPhases::Aqua]? rates[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
         const auto& liquid = pu.phase_used[BlackoilPhases::Liquid]? rates[pu.phase_pos[BlackoilPhases::Liquid]]:0.0;
@@ -82,7 +82,7 @@ public:
     }
 
     double getExplicitGFR(const int table_id, const size_t well_index) const {
-        const auto& rates = well_state_.well(well_index).surface_rates;
+        const auto& rates = well_state_.well(well_index).prev_surface_rates;
         const auto& pu = well_state_.phaseUsage();
         const auto& aqua = pu.phase_used[BlackoilPhases::Aqua]? rates[pu.phase_pos[BlackoilPhases::Aqua]]:0.0;
         const auto& liquid = pu.phase_used[BlackoilPhases::Liquid]? rates[pu.phase_pos[BlackoilPhases::Liquid]]:0.0;

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -347,6 +347,17 @@ void WellInterfaceGeneric::setVFPProperties(const VFPProperties* vfp_properties_
     vfp_properties_ = vfp_properties_arg;
 }
 
+void WellInterfaceGeneric::setPrevSurfaceRates(WellState& well_state, 
+                                               const WellState& prev_well_state) const
+    {
+        auto& ws = well_state.well(this->index_of_well_);
+        if (!this->changedToOpenThisStep()){
+            ws.prev_surface_rates = prev_well_state.well(this->index_of_well_).surface_rates;
+        } else {
+            ws.prev_surface_rates = ws.surface_rates;
+        }
+    }
+
 void WellInterfaceGeneric::setGuideRate(const GuideRate* guide_rate_arg)
 {
     guide_rate_ = guide_rate_arg;
@@ -471,7 +482,7 @@ bool WellInterfaceGeneric::isOperableAndSolvable() const
 bool WellInterfaceGeneric::useVfpExplicit() const
 {
     const auto& wvfpexp = well_ecl_.getWVFPEXP();
-    return ((wvfpexp.explicit_lookup() && !changedToOpenThisStep())|| operability_status_.use_vfpexplicit);
+    return (wvfpexp.explicit_lookup() || operability_status_.use_vfpexplicit);
 }
 
 bool WellInterfaceGeneric::thpLimitViolatedButNotSwitched() const

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -94,6 +94,7 @@ public:
     void closeCompletions(const WellTestState& wellTestState);
 
     void setVFPProperties(const VFPProperties* vfp_properties_arg);
+    void setPrevSurfaceRates(WellState& well_state, const WellState& prev_well_state) const;
     void setGuideRate(const GuideRate* guide_rate_arg);
     void setWellEfficiencyFactor(const double efficiency_factor);
     void setRepRadiusPerfLength();

regressionTests.cmake

@@ -1191,14 +1191,16 @@ add_test_compareECLFiles(CASENAME 01_wgrupcon
                          SIMULATOR flow
                          ABS_TOL ${abs_tol}
                          REL_TOL ${rel_tol}
-                         DIR wgrupcon)
+                         DIR wgrupcon
+                         TEST_ARGS --enable-tuning=true)
 
 add_test_compareECLFiles(CASENAME 02_wgrupcon
                          FILENAME 02-WGRUPCON
                          SIMULATOR flow
                          ABS_TOL ${abs_tol}
                          REL_TOL ${rel_tol}
-                         DIR wgrupcon)
+                         DIR wgrupcon
+                         TEST_ARGS --enable-tuning=true)
 add_test_compareECLFiles(CASENAME winjmult_stdw
                          FILENAME WINJMULT_STDW
                          SIMULATOR flow