adding function checkConvergenceControlEq to refactor getWellConvergence

only functional change is that rate control or BHP control has different
tolerance for MSW now.

opm/simulators/wells/MultisegmentWell.hpp

@@ -379,7 +379,13 @@ namespace Opm
         void detectOscillations(const std::vector<double>& measure_history,
                                 const int it, bool& oscillate, bool& stagnate) const;
 
-        double getResidualMeasureValue(const std::vector<double>& residuals) const;
+        double getResidualMeasureValue(const std::vector<double>& residuals,
+                                       DeferredLogger& deferred_logger) const;
+
+        double getControlTolerance(DeferredLogger& deferred_logger) const;
+
+        void checkConvergenceControlEq(ConvergenceReport& report,
+                                       DeferredLogger& deferred_logger) const;
 
     };
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -429,28 +429,9 @@ namespace Opm
             }
         }
 
-        using CR = ConvergenceReport;
-        CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
-        switch(well_controls_get_current_type(well_controls_)) {
-            case THP:
-                ctrltype = CR::WellFailure::Type::ControlTHP;
-                break;
-            case BHP:
-                ctrltype = CR::WellFailure::Type::ControlBHP;
-                break;
-            case RESERVOIR_RATE:
-            case SURFACE_RATE:
-                ctrltype = CR::WellFailure::Type::ControlRate;
-                break;
-            default:
-                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
-        }
-        assert(ctrltype != CR::WellFailure::Type::Invalid);
-
         std::vector<double> maximum_residual(numWellEq, 0.0);
 
         ConvergenceReport report;
-        const int dummy_component = -1;
         // TODO: the following is a little complicated, maybe can be simplified in some way?
         for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
             for (int seg = 0; seg < numberOfSegments(); ++seg) {
@@ -460,18 +441,8 @@ namespace Opm
                         maximum_residual[eq_idx] = flux_residual;
                     }
                 } else { // pressure or control equation
-                    if (seg == 0) {
-                        // Control equation
-                        const double control_residual = abs_residual[seg][eq_idx];
-                        if (std::isnan(control_residual)) {
-                            report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
-                        } else if (control_residual > param_.max_residual_allowed_) {
-                            report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
-                            // TODO: we should distinguish the flux residual or pressure residual here
-                        } else if (control_residual > param_.tolerance_wells_) {
-                            report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
-                        }
-                    } else {
+                    // for the top segment (seg == 0), it is control equation, will be checked later separately
+                    if (seg > 0) {
                         // Pressure equation
                         const double pressure_residual = abs_residual[seg][eq_idx];
                         if (pressure_residual > maximum_residual[eq_idx]) {
@@ -482,6 +453,7 @@ namespace Opm
             }
         }
 
+        using CR = ConvergenceReport;
         for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
             if (eq_idx < num_components_) { // phase or component mass equations
                 const double flux_residual = maximum_residual[eq_idx];
@@ -495,6 +467,7 @@ namespace Opm
                 }
             } else { // pressure equation
                 const double pressure_residual = maximum_residual[eq_idx];
+                const int dummy_component = -1;
                 if (std::isnan(pressure_residual)) {
                     report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, name()});
                 } else if (std::isinf(pressure_residual)) {
@@ -505,6 +478,8 @@ namespace Opm
             }
         }
 
+        checkConvergenceControlEq(report, deferred_logger);
+
         return report;
     }
 
@@ -1825,7 +1800,7 @@ namespace Opm
             }
 
             residual_history.push_back(getWellResiduals(B_avg));
-            measure_history.push_back(getResidualMeasureValue(residual_history[it]));
+            measure_history.push_back(getResidualMeasureValue(residual_history[it], deferred_logger) );
 
             bool is_oscillate = false;
             bool is_stagnate = false;
@@ -2248,7 +2223,8 @@ namespace Opm
     template<typename TypeTag>
     double
     MultisegmentWell<TypeTag>::
-    getResidualMeasureValue(const std::vector<double>& residuals) const
+    getResidualMeasureValue(const std::vector<double>& residuals,
+                            DeferredLogger& deferred_logger) const
     {
         assert(int(residuals.size()) == numWellEq + 1);
 
@@ -2262,9 +2238,11 @@ namespace Opm
             }
         }
 
-        const double pressure_tolerance = param_.tolerance_pressure_ms_wells_;
-        if (residuals[SPres] > pressure_tolerance) {
-            sum += residuals[SPres] / pressure_tolerance;
+        const double control_tolerance = getControlTolerance(deferred_logger);
+
+        // const double pressure_tolerance = param_.tolerance_pressure_ms_wells_;
+        if (residuals[SPres] > control_tolerance) {
+            sum += residuals[SPres] / control_tolerance;
             ++count;
         }
 
@@ -2275,4 +2253,69 @@ namespace Opm
         return sum;
     }
 
+
+
+
+
+    template<typename TypeTag>
+    double
+    MultisegmentWell<TypeTag>::
+    getControlTolerance(DeferredLogger& deferred_logger) const
+    {
+        double control_tolerance = 0.;
+        switch(well_controls_get_current_type(well_controls_) ) {
+            case BHP:
+            case THP:
+                control_tolerance = param_.tolerance_wells_;
+                break;
+            case RESERVOIR_RATE:
+            case SURFACE_RATE:
+                control_tolerance = param_.tolerance_pressure_ms_wells_;
+                break;
+            default:
+                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
+        }
+
+        return control_tolerance;
+    }
+
+
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    checkConvergenceControlEq(ConvergenceReport& report,
+                              DeferredLogger& deferred_logger) const
+    {
+        using CR = ConvergenceReport;
+        CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
+        switch(well_controls_get_current_type(well_controls_)) {
+            case THP:
+                ctrltype = CR::WellFailure::Type::ControlTHP;
+                break;
+            case BHP:
+                ctrltype = CR::WellFailure::Type::ControlBHP;
+                break;
+            case RESERVOIR_RATE:
+            case SURFACE_RATE:
+                ctrltype = CR::WellFailure::Type::ControlRate;
+                break;
+            default:
+                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
+        }
+        assert(ctrltype != CR::WellFailure::Type::Invalid);
+
+        const double control_residual = std::abs(resWell_[0][SPres]);
+        const double control_tolerance = getControlTolerance(deferred_logger);
+
+        const int dummy_component = -1;
+        if (std::isnan(control_residual)) {
+            report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
+        } else if (control_residual > param_.max_residual_allowed_) {
+            report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
+        } else if (control_residual > control_tolerance) {
+            report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
+        }
+    }
 }

opm/simulators/wells/StandardWell.hpp

@@ -439,6 +439,10 @@ namespace Opm
                                          Opm::DeferredLogger& deferred_logger) override;
 
         virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
+
+        void checkConvergenceControlEq(ConvergenceReport& report,
+                                       DeferredLogger& deferred_logger) const;
+
     };
 
 }

opm/simulators/wells/StandardWellV.hpp

@@ -465,6 +465,9 @@ namespace Opm
                                                Opm::DeferredLogger& deferred_logger);
 
         virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
+
+        void checkConvergenceControlEq(ConvergenceReport& report,
+                                       DeferredLogger& deferred_logger) const;
     };
 
 }

opm/simulators/wells/StandardWellV_impl.hpp

@@ -2029,40 +2029,12 @@ namespace Opm
             }
         }
 
-        // processing the residual of the well control equation
-        const double well_control_residual = res[numWellEq_ - 1];
-        // TODO: we should have better way to specify the control equation tolerance
-        double control_tolerance = 0.;
-        switch(well_controls_get_current_type(well_controls_)) {
-            case THP:
-                type = CR::WellFailure::Type::ControlTHP;
-                control_tolerance = 1.e4; // 0.1 bar
-                break;
-            case BHP:  // pressure type of control
-                type = CR::WellFailure::Type::ControlBHP;
-                control_tolerance = 1.e3; // 0.01 bar
-                break;
-            case RESERVOIR_RATE:
-            case SURFACE_RATE:
-                type = CR::WellFailure::Type::ControlRate;
-                control_tolerance = 1.e-4; // smaller tolerance for rate control
-                break;
-            default:
-                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " +  name(), deferred_logger);
-        }
-
-        const int dummy_component = -1;
-        if (std::isnan(well_control_residual)) {
-            report.setWellFailed({type, CR::Severity::NotANumber, dummy_component, name()});
-        } else if (well_control_residual > maxResidualAllowed * 10.) {
-            report.setWellFailed({type, CR::Severity::TooLarge, dummy_component, name()});
-        } else if ( well_control_residual > control_tolerance) {
-            report.setWellFailed({type, CR::Severity::Normal, dummy_component, name()});
-        }
+        checkConvergenceControlEq(report, deferred_logger);
 
         if (this->has_polymermw && well_type_ == INJECTOR) {
             //  checking the convergence of the perforation rates
             const double wat_vel_tol = 1.e-8;
+            const int dummy_component = -1;
             const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
             for (int perf = 0; perf < number_of_perforations_; ++perf) {
                 const double wat_vel_residual = res[Bhp + 1 + perf];
@@ -3112,4 +3084,49 @@ namespace Opm
             duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
         }
     }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWellV<TypeTag>::
+    checkConvergenceControlEq(ConvergenceReport& report,
+                              DeferredLogger& deferred_logger) const
+    {
+        double control_tolerance = 0.;
+        using CR = ConvergenceReport;
+        CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
+        switch(well_controls_get_current_type(well_controls_)) {
+            case THP:
+                ctrltype = CR::WellFailure::Type::ControlTHP;
+                control_tolerance = 1.e4; // 0.1 bar
+                break;
+            case BHP:  // pressure type of control
+                ctrltype = CR::WellFailure::Type::ControlBHP;
+                control_tolerance = 1.e3; // 0.01 bar
+                break;
+            case RESERVOIR_RATE:
+            case SURFACE_RATE:
+                ctrltype = CR::WellFailure::Type::ControlRate;
+                control_tolerance = 1.e-4; // smaller tolerance for rate control
+                break;
+            default:
+                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
+        }
+        assert(ctrltype != CR::WellFailure::Type::Invalid);
+
+        const double well_control_residual = std::abs(resWell_[0][Bhp]);
+        const int dummy_component = -1;
+        const double max_residual_allowed = param_.max_residual_allowed_;
+        if (std::isnan(well_control_residual)) {
+            report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
+        } else if (well_control_residual > max_residual_allowed * 10.) {
+            report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
+        } else if ( well_control_residual > control_tolerance) {
+            report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
+        }
+
+    }
 }

opm/simulators/wells/StandardWell_impl.hpp

@@ -1957,37 +1957,7 @@ namespace Opm
             }
         }
 
-        // processing the residual of the well control equation
-        const double well_control_residual = res[numWellEq - 1];
-        // TODO: we should have better way to specify the control equation tolerance
-        double control_tolerance = 0.;
-        switch(well_controls_get_current_type(well_controls_)) {
-            case THP:
-                type = CR::WellFailure::Type::ControlTHP;
-                control_tolerance = 1.e4; // 0.1 bar
-                break;
-            case BHP:  // pressure type of control
-                type = CR::WellFailure::Type::ControlBHP;
-                control_tolerance = 1.e3; // 0.01 bar
-                break;
-            case RESERVOIR_RATE:
-            case SURFACE_RATE:
-                type = CR::WellFailure::Type::ControlRate;
-                control_tolerance = 1.e-4; // smaller tolerance for rate control
-                break;
-            default:
-                // If this happens, control_tolerance will remain zero and convergence will fail.
-                deferred_logger.bug("Unknown well control type for well " +  name());
-        }
-
-        const int dummy_component = -1;
-        if (std::isnan(well_control_residual)) {
-            report.setWellFailed({type, CR::Severity::NotANumber, dummy_component, name()});
-        } else if (well_control_residual > maxResidualAllowed * 10.) {
-            report.setWellFailed({type, CR::Severity::TooLarge, dummy_component, name()});
-        } else if ( well_control_residual > control_tolerance) {
-            report.setWellFailed({type, CR::Severity::Normal, dummy_component, name()});
-        }
+        checkConvergenceControlEq(report, deferred_logger);
 
         return report;
     }
@@ -2853,4 +2823,48 @@ namespace Opm
     updateWaterThroughput(const double dt OPM_UNUSED, WellState& well_state OPM_UNUSED) const
     {
     }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    checkConvergenceControlEq(ConvergenceReport& report,
+                              DeferredLogger& deferred_logger) const
+    {
+        double control_tolerance = 0.;
+        using CR = ConvergenceReport;
+        CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
+        switch(well_controls_get_current_type(well_controls_)) {
+            case THP:
+                ctrltype = CR::WellFailure::Type::ControlTHP;
+                control_tolerance = 1.e4; // 0.1 bar
+                break;
+            case BHP:  // pressure type of control
+                ctrltype = CR::WellFailure::Type::ControlBHP;
+                control_tolerance = 1.e3; // 0.01 bar
+                break;
+            case RESERVOIR_RATE:
+            case SURFACE_RATE:
+                ctrltype = CR::WellFailure::Type::ControlRate;
+                control_tolerance = 1.e-4; // smaller tolerance for rate control
+                break;
+            default:
+                OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << name(), deferred_logger);
+        }
+        assert(ctrltype != CR::WellFailure::Type::Invalid);
+
+        const double well_control_residual = std::abs(resWell_[0][numWellEq -1]);
+        const int dummy_component = -1;
+        const double max_residual_allowed = param_.max_residual_allowed_;
+        if (std::isnan(well_control_residual)) {
+            report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()});
+        } else if (well_control_residual > max_residual_allowed * 10.) {
+            report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()});
+        } else if ( well_control_residual > control_tolerance) {
+            report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()});
+        }
+    }
 }