Merge pull request #5822 from totto82/fixNETV4

Fixes related to zero guiderates from NETV

opm/simulators/wells/FractionCalculator.cpp

@@ -77,18 +77,25 @@ Scalar FractionCalculator<Scalar>::
 localFraction(const std::string& name,
               const std::string& always_included_child)
 {
-    const Scalar my_guide_rate = guideRate(name, always_included_child);
+    bool always_use_potentials = false;
-    const Group& parent_group = schedule_.getGroup(parent(name), report_step_);
+    const Scalar my_guide_rate = guideRate(name, always_included_child, always_use_potentials);
-    const Scalar total_guide_rate = guideRateSum(parent_group, always_included_child);
 
-    // the total guide gate is the same as my_guide rate
+    const Group& parent_group = schedule_.getGroup(parent(name), report_step_);
-    // the well/group is probably on its own, i.e. return 1
+    const auto [total_guide_rate, num_active_groups] = guideRateSum(parent_group, always_included_child, always_use_potentials);
-    // even is its guide_rate is zero
+
-    const Scalar guide_rate_epsilon = 1e-12;
+    // the group/well "name" is the only active group/well we therefore return 1 as the fraction
-    if ( std::abs(my_guide_rate - total_guide_rate) < guide_rate_epsilon )
+    // even though my_guide_rate may be zero
+    if (num_active_groups == 1)
         return 1.0;
 
-    assert(total_guide_rate > my_guide_rate);
+    if (total_guide_rate == 0 ) {
+        // if the total guide rate is zero (for instance due to netv = 0) we use the potentials
+        // to distribute the group rate
+        always_use_potentials = true;
+        const Scalar my_pot = guideRate(name, always_included_child, always_use_potentials);
+        const Scalar my_total_pot = guideRateSum(parent_group, always_included_child, always_use_potentials).first;
+        return my_pot / my_total_pot;
+    }
     return my_guide_rate / total_guide_rate;
 }
 
@@ -104,11 +111,13 @@ parent(const std::string& name)
 }
 
 template<class Scalar>
-Scalar FractionCalculator<Scalar>::
+std::pair<Scalar, int> FractionCalculator<Scalar>::
 guideRateSum(const Group& group,
-             const std::string& always_included_child)
+             const std::string& always_included_child,
+             const bool always_use_potentials)
 {
     Scalar total_guide_rate = 0.0;
+    int number_of_included_well_or_groups = 0;
     for (const std::string& child_group : group.groups()) {
         bool included = (child_group == always_included_child);
         if (is_producer_) {
@@ -122,7 +131,10 @@ guideRateSum(const Group& group,
                         (ctrl == Group::InjectionCMode::NONE);
         }
         if (included) {
-            total_guide_rate += guideRate(child_group, always_included_child);
+            if (groupControlledWells(child_group, always_included_child) > 0) {
+                number_of_included_well_or_groups++;
+                total_guide_rate += guideRate(child_group, always_included_child, always_use_potentials);
+            }
         }
     }
 
@@ -133,34 +145,35 @@ guideRateSum(const Group& group,
         } else {
             included |= well_state_.isInjectionGrup(child_well);
         }
-
         if (included) {
-            total_guide_rate += guideRate(child_well, always_included_child);
+            number_of_included_well_or_groups++;
+            total_guide_rate += guideRate(child_well, always_included_child, always_use_potentials);
         }
     }
-    return total_guide_rate;
+    return {total_guide_rate, number_of_included_well_or_groups};
 }
 
 template<class Scalar>
 Scalar FractionCalculator<Scalar>::
 guideRate(const std::string& name,
-          const std::string& always_included_child)
+          const std::string& always_included_child,
+          const bool always_use_potentials)
 {
     if (schedule_.hasWell(name, report_step_)) {
         return WellGroupHelpers<Scalar>::getGuideRate(name, schedule_, well_state_, group_state_,
                                                       report_step_, guide_rate_, target_, pu_);
     } else {
         if (groupControlledWells(name, always_included_child) > 0) {
-            if (is_producer_ && guide_rate_->has(name)) {
+            if (is_producer_ && guide_rate_->has(name) && !always_use_potentials) {
                 return guide_rate_->get(name, target_, getGroupRateVector(name));
-            } else if (!is_producer_ && guide_rate_->has(name, injection_phase_)) {
+            } else if (!is_producer_ && guide_rate_->has(name, injection_phase_) && !always_use_potentials) {
                 return guide_rate_->get(name, injection_phase_);
             } else {
                 // We are a group, with default guide rate.
                 // Compute guide rate by accumulating our children's guide rates.
                 const Group& group = schedule_.getGroup(name, report_step_);
                 const double eff = group.getGroupEfficiencyFactor();
-                return eff * guideRateSum(group, always_included_child);
+                return eff * guideRateSum(group, always_included_child, always_use_potentials).first;
             }
         } else {
             // No group-controlled subordinate wells.

opm/simulators/wells/FractionCalculator.hpp

@@ -56,10 +56,15 @@ public:
 
 private:
     std::string parent(const std::string& name);
-    Scalar guideRateSum(const Group& group,
+
-                        const std::string& always_included_child);
+    // returns the sum of the guiderates of the given group
+    // and the number of sub-groups/wells that contributed to the sum
+    std::pair<Scalar,int> guideRateSum(const Group& group,
+                        const std::string& always_included_child,
+                        const bool always_use_potentials);
     Scalar guideRate(const std::string& name,
-                     const std::string& always_included_child);
+                     const std::string& always_included_child,
+                     const bool always_use_potentials);
     int groupControlledWells(const std::string& group_name,
                              const std::string& always_included_child);
     GuideRate::RateVector getGroupRateVector(const std::string& group_name);

opm/simulators/wells/StandardWell_impl.hpp

@@ -1758,7 +1758,10 @@ namespace Opm
         }
 
         bool converged_implicit = false;
-        if (this->param_.local_well_solver_control_switching_) {
+        // for newly opened wells we dont compute the potentials implicit
+        // group controlled wells with defaulted guiderates will have zero targets as
+        // the potentials are used to compute the well fractions.
+        if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) {
             converged_implicit = computeWellPotentialsImplicit(simulator, well_potentials, deferred_logger);
         }
         if (!converged_implicit) {

opm/simulators/wells/WellGroupHelpers.cpp

@@ -300,14 +300,14 @@ updateGuideRatesForInjectionGroups(const Group& group,
         if(!group.hasInjectionControl(phase))
             continue;
 
-        Scalar guideRateValue = 0.0;
+        std::optional<Scalar> guideRateValue;
         const auto& controls = group.injectionControls(phase, summaryState);
         switch (controls.guide_rate_def){
         case Group::GuideRateInjTarget::RATE:
             break;
         case Group::GuideRateInjTarget::VOID:
         {
-            guideRateValue = group_state.injection_vrep_rate(group.name());
+            guideRateValue = std::max(Scalar(0.0), group_state.injection_vrep_rate(group.name()));
             break;
         }
         case Group::GuideRateInjTarget::NETV:
@@ -315,17 +315,17 @@ updateGuideRatesForInjectionGroups(const Group& group,
             guideRateValue = group_state.injection_vrep_rate(group.name());
             const std::vector<Scalar>& injRES = group_state.injection_reservoir_rates(group.name());
             if (phase != Phase::OIL && pu.phase_used[BlackoilPhases::Liquid])
-                guideRateValue -= injRES[pu.phase_pos[BlackoilPhases::Liquid]];
+                guideRateValue = *guideRateValue - injRES[pu.phase_pos[BlackoilPhases::Liquid]];
             if (phase != Phase::GAS && pu.phase_used[BlackoilPhases::Vapour])
-                guideRateValue -= injRES[pu.phase_pos[BlackoilPhases::Vapour]];
+                guideRateValue = *guideRateValue - injRES[pu.phase_pos[BlackoilPhases::Vapour]];
             if (phase != Phase::WATER && pu.phase_used[BlackoilPhases::Aqua])
-                guideRateValue -= injRES[pu.phase_pos[BlackoilPhases::Aqua]];
+                guideRateValue = *guideRateValue - injRES[pu.phase_pos[BlackoilPhases::Aqua]];
+
+            guideRateValue = std::max(Scalar(0.0), *guideRateValue);
             break;
         }
         case Group::GuideRateInjTarget::RESV:
             OPM_DEFLOG_THROW(std::runtime_error, "GUIDE PHASE RESV not implemented. Group " + group.name(), deferred_logger);
-        case Group::GuideRateInjTarget::POTN:
-            break;
         case Group::GuideRateInjTarget::NO_GUIDE_RATE:
             break;
         default:
@@ -335,7 +335,9 @@ updateGuideRatesForInjectionGroups(const Group& group,
         }
 
         const UnitSystem& unit_system = schedule.getUnits();
-        guideRateValue = unit_system.from_si(UnitSystem::measure::rate, guideRateValue);
+        if (guideRateValue) {
+            guideRateValue = unit_system.from_si(UnitSystem::measure::rate, *guideRateValue);
+        }
         guideRate->compute(group.name(), phase, reportStepIdx, guideRateValue);
     }
 }

regressionTests.cmake

@@ -954,7 +954,8 @@ add_test_compareECLFiles(CASENAME 9_3d_grpctl_stw_model2
                          SIMULATOR flow
                          ABS_TOL ${abs_tol}
                          REL_TOL ${rel_tol}
-                         DIR model2)
+                         DIR model2
+                         TEST_ARGS --time-step-after-event-in-days=1)
 
 add_test_compareECLFiles(CASENAME 9_3d_grpctl_msw_model2
                          FILENAME 9_3D_GINJ_GAS_MAX_EXPORT_MSW