Expose hardcoded well/Newton tuning constants; remove the dead tolerance-well-control knob

flowexperimental/flowexp.hpp

@@ -160,7 +160,6 @@ class FlowExpProblem : public FlowProblemBlackoil<TypeTag> //, public FvBaseProb
         Parameters::Hide<Parameters::ToleranceCnv<Scalar>>();
         Parameters::Hide<Parameters::ToleranceCnvRelaxed<Scalar>>();
         Parameters::Hide<Parameters::ToleranceWells<Scalar>>();
-        Parameters::Hide<Parameters::ToleranceWellControl<Scalar>>();
         Parameters::Hide<Parameters::MaxWelleqIter>();
         Parameters::Hide<Parameters::UseMultisegmentWell>();
         Parameters::Hide<Parameters::TolerancePressureMsWells<Scalar>>();

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -56,10 +56,14 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     tolerance_max_drs_ = Parameters::Get<Parameters::ToleranceMaxDrs<Scalar>>();
     tolerance_max_drv_ = Parameters::Get<Parameters::ToleranceMaxDrv<Scalar>>();
     tolerance_wells_ = Parameters::Get<Parameters::ToleranceWells<Scalar>>();
-    tolerance_well_control_ = Parameters::Get<Parameters::ToleranceWellControl<Scalar>>();
+    tolerance_wells_stopped_factor_ = Parameters::Get<Parameters::ToleranceWellsStoppedFactor<Scalar>>();
+    tolerance_wells_dynamic_thp_factor_ = Parameters::Get<Parameters::ToleranceWellsDynamicThpFactor<Scalar>>();
     max_welleq_iter_ = Parameters::Get<Parameters::MaxWelleqIter>();
     use_multisegment_well_ = Parameters::Get<Parameters::UseMultisegmentWell>();
     tolerance_pressure_ms_wells_ = Parameters::Get<Parameters::TolerancePressureMsWells<Scalar>>();
+    tolerance_well_bhp_eq_ = Parameters::Get<Parameters::ToleranceWellBhpEq<Scalar>>();
+    tolerance_well_thp_eq_ = Parameters::Get<Parameters::ToleranceWellThpEq<Scalar>>();
+    tolerance_well_grup_eq_ = Parameters::Get<Parameters::ToleranceWellGrupEq<Scalar>>();
     relaxed_tolerance_flow_well_ = Parameters::Get<Parameters::RelaxedWellFlowTol<Scalar>>();
     relaxed_tolerance_pressure_ms_well_ = Parameters::Get<Parameters::RelaxedPressureTolMsw<Scalar>>();
     max_pressure_change_ms_wells_ = Parameters::Get<Parameters::MaxPressureChangeMsWells<Scalar>>();
@@ -70,6 +74,8 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     max_niter_inner_well_iter_ = Parameters::Get<Parameters::MaxNewtonIterationsWithInnerWellIterations>();
     shut_unsolvable_wells_ = Parameters::Get<Parameters::ShutUnsolvableWells>();
     max_inner_iter_wells_ = Parameters::Get<Parameters::MaxInnerIterWells>();
+    min_iter_after_switch_wells_ = Parameters::Get<Parameters::MinIterAfterSwitchWells>();
+    min_iter_after_switch_ms_wells_ = Parameters::Get<Parameters::MinIterAfterSwitchMsWells>();
     max_well_status_switch_inner_iter_ = Parameters::Get<Parameters::MaxWellStatusSwitchInInnerIterWells>();
     max_well_status_switch_ = Parameters::Get<Parameters::MaxWellStatusSwitchForWells>();
     maxSinglePrecisionTimeStep_ = Parameters::Get<Parameters::MaxSinglePrecisionDays<Scalar>>() * 24 * 60 * 60;
@@ -185,8 +191,12 @@ void BlackoilModelParameters<Scalar>::registerParameters()
          "of residual tolerances. Use with care!");
     Parameters::Register<Parameters::ToleranceWells<Scalar>>
         ("Well convergence tolerance");
-    Parameters::Register<Parameters::ToleranceWellControl<Scalar>>
-        ("Tolerance for the well control equations");
+    Parameters::Register<Parameters::ToleranceWellsStoppedFactor<Scalar>>
+        ("Multiplier applied to the well convergence tolerance for stopped wells "
+         "and wells under a zero rate target");
+    Parameters::Register<Parameters::ToleranceWellsDynamicThpFactor<Scalar>>
+        ("Multiplier applied to the well convergence tolerance for wells on a "
+         "dynamic THP limit, to help network convergence");
     Parameters::Register<Parameters::MaxWelleqIter>
         ("Maximum number of iterations to determine solution the well equations");
     Parameters::Register<Parameters::UseMultisegmentWell>
@@ -198,6 +208,12 @@ void BlackoilModelParameters<Scalar>::registerParameters()
          "or 'per-connection' (one linear tubing, a segment per connection)");
     Parameters::Register<Parameters::TolerancePressureMsWells<Scalar>>
         ("Tolerance for the pressure equations for multi-segment wells");
+    Parameters::Register<Parameters::ToleranceWellBhpEq<Scalar>>
+        ("Tolerance for the BHP control equation residual of standard wells [Pa]");
+    Parameters::Register<Parameters::ToleranceWellThpEq<Scalar>>
+        ("Tolerance for the THP control equation residual of standard wells [Pa]");
+    Parameters::Register<Parameters::ToleranceWellGrupEq<Scalar>>
+        ("Tolerance for the group-control equation residual of standard wells [m3/s]");
     Parameters::Register<Parameters::RelaxedWellFlowTol<Scalar>>
         ("Relaxed tolerance for the well flow residual");
     Parameters::Register<Parameters::RelaxedPressureTolMsw<Scalar>>
@@ -217,6 +233,12 @@ void BlackoilModelParameters<Scalar>::registerParameters()
         ("Shut unsolvable wells");
     Parameters::Register<Parameters::MaxInnerIterWells>
         ("Maximum number of inner iterations for standard wells");
+    Parameters::Register<Parameters::MinIterAfterSwitchWells>
+        ("Inner iterations forced between two control switches of a standard well. "
+         "Values <= 1 let the inner solve cycle between controls.");
+    Parameters::Register<Parameters::MinIterAfterSwitchMsWells>
+        ("Inner iterations forced between two control switches of a multi-segment well. "
+         "Values <= 1 let the inner solve cycle between controls.");
     Parameters::Register<Parameters::MaxWellStatusSwitchInInnerIterWells>
         ("Maximum number of status switching (stop<->open) for a well during inner iterations.");
     Parameters::Register<Parameters::MaxWellStatusSwitchForWells>
@@ -230,7 +252,9 @@ void BlackoilModelParameters<Scalar>::registerParameters()
          "arithmetic can be used solving for the linear systems of equations");
     Parameters::Register<Parameters::MinStrictCnvIter>
         ("Minimum number of Newton iterations before relaxed tolerances "
-         "can be used for the CNV convergence criterion");
+         "can be used for the CNV convergence criterion. "
+         "Default -1 means that the relaxed tolerance is used when maximum "
+         "number of Newton iterations are reached.");
     Parameters::Register<Parameters::MinStrictMbIter>
         ("Minimum number of Newton iterations before relaxed tolerances "
          "can be used for the MB convergence criterion. "

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -91,7 +91,10 @@ template<class Scalar>
 struct ToleranceWells { static constexpr Scalar value = 1e-4; };
 
 template<class Scalar>
-struct ToleranceWellControl { static constexpr Scalar value = 1e-7; };
+struct ToleranceWellsStoppedFactor { static constexpr Scalar value = 1e-4; };
+
+template<class Scalar>
+struct ToleranceWellsDynamicThpFactor { static constexpr Scalar value = 1e-1; };
 
 struct MaxWelleqIter { static constexpr int value = 30; };
 
@@ -116,13 +119,23 @@ struct ConvertToMultisegmentWell { static constexpr auto value = "none"; };
 
 template<class Scalar>
 struct TolerancePressureMsWells { static constexpr Scalar value = 0.01*1e5; };
+/// Control-equation tolerances for standard wells. Defaults reproduce the values that
+/// used to be hardcoded at the call site; MSW has had these as parameters all along.
+template<class Scalar>
+struct ToleranceWellBhpEq { static constexpr Scalar value = 1.0e3; };
+template<class Scalar>
+struct ToleranceWellThpEq { static constexpr Scalar value = 1.0e4; };
+template<class Scalar>
+struct ToleranceWellGrupEq { static constexpr Scalar value = 1.0e-6; };
 
 template<class Scalar>
 struct MaxPressureChangeMsWells { static constexpr Scalar value = 10*1e5; };
 
 struct MaxNewtonIterationsWithInnerWellIterations { static constexpr int value = 99; };
 struct MaxInnerIterMsWells { static constexpr int value = 100; };
 struct MaxInnerIterWells { static constexpr int value = 50; };
+struct MinIterAfterSwitchWells { static constexpr int value = 4; };
+struct MinIterAfterSwitchMsWells { static constexpr int value = 3; };
 struct MaxWellStatusSwitchInInnerIterWells { static constexpr int value = 99; };
 struct MaxWellStatusSwitchForWells { static constexpr int value = 99; };
 struct ShutUnsolvableWells { static constexpr bool value = true; };
@@ -247,11 +260,16 @@ struct BlackoilModelParameters
     Scalar tolerance_max_drv_;
     /// Well convergence tolerance.
     Scalar tolerance_wells_;
-    /// Tolerance for the well control equations
-    //  TODO: it might need to distinguish between rate control and pressure control later
-    Scalar tolerance_well_control_;
+    /// Multipliers applied to tolerance_wells_ for stopped/zero-rate wells and for
+    /// wells on a dynamic THP limit (the latter to help network convergence).
+    Scalar tolerance_wells_stopped_factor_;
+    Scalar tolerance_wells_dynamic_thp_factor_;
     /// Tolerance for the pressure equations for multisegment wells
     Scalar tolerance_pressure_ms_wells_;
+    /// standard-well control-equation tolerances (BHP/THP in Pa, GRUP in m3/s)
+    Scalar tolerance_well_bhp_eq_;
+    Scalar tolerance_well_thp_eq_;
+    Scalar tolerance_well_grup_eq_;
     /// Relaxed tolerance for for the well flow residual
     Scalar relaxed_tolerance_flow_well_;
 
@@ -282,6 +300,11 @@ struct BlackoilModelParameters
     /// Maximum inner iteration number for standard wells
     int max_inner_iter_wells_;
 
+    /// Inner iterations forced between two control switches of the same well.
+    /// Values <= 1 are known to let the inner solve cycle between controls.
+    int min_iter_after_switch_wells_;
+    int min_iter_after_switch_ms_wells_;
+
     /// Maximum iteration number of the well equation solution
     int max_welleq_iter_;
 

opm/simulators/flow/NonlinearSolver.cpp

@@ -146,6 +146,8 @@ NonlinearSolverParameters()
 
     // overload with given parameters
     relaxMax_ = Parameters::Get<Parameters::NewtonMaxRelax<Scalar>>();
+    relaxIncrement_ = Parameters::Get<Parameters::NewtonRelaxIncrement<Scalar>>();
+    relaxRelTol_ = Parameters::Get<Parameters::NewtonRelaxRelTol<Scalar>>();
 
     const auto& relaxationTypeString = Parameters::Get<Parameters::NewtonRelaxationType>();
     if (relaxationTypeString == "dampen") {
@@ -175,6 +177,12 @@ registerParameters()
 {
     Parameters::Register<Parameters::NewtonMaxRelax<Scalar>>
         ("The maximum relaxation factor of a Newton iteration");
+    Parameters::Register<Parameters::NewtonRelaxIncrement<Scalar>>
+        ("Amount by which the Newton relaxation factor is reduced each time an "
+         "oscillation is detected");
+    Parameters::Register<Parameters::NewtonRelaxRelTol<Scalar>>
+        ("Relative tolerance used to classify a residual history as oscillating "
+         "rather than stagnating");
     Parameters::Register<Parameters::NewtonRelaxationType>
         ("The type of relaxation used by Newton method. Valid options are: dampen or sor");
 }

opm/simulators/flow/NonlinearSolver.hpp

@@ -45,6 +45,12 @@ namespace Opm::Parameters {
 template<class Scalar>
 struct NewtonMaxRelax { static constexpr Scalar value = 0.5; };
 
+template<class Scalar>
+struct NewtonRelaxIncrement { static constexpr Scalar value = 0.1; };
+
+template<class Scalar>
+struct NewtonRelaxRelTol { static constexpr Scalar value = 0.2; };
+
 struct NewtonRelaxationType { static constexpr auto value = "dampen"; };
 
 } // namespace Opm::Parameters

opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp

@@ -810,8 +810,8 @@ getReservoirConvergence(const double reportTime,
         || relax_iter_cnv
         || relax_dsol_cnv;
 
-    // Ensure that CNV convergence criteria is met when max.
-    // solution change tolerances have been fulfilled
+    // The solution-change criterion overrides CNV entirely rather than merely relaxing it:
+    // 1e20 is never exceeded, so a converged solution update accepts any CNV.
     Scalar tolerance_cnv_relaxed = relax_dsol_cnv ? 1e20 : this->param_.tolerance_cnv_relaxed_;
 
     const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : this->param_.tolerance_cnv_;

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1705,7 +1705,7 @@ namespace Opm
         // Always take a few (more than one) iterations after a switch before allowing a new switch
         // The optimal number here is subject to further investigation, but it has been observerved
         // that unless this number is >1, we may get stuck in a cycle
-        const int min_its_after_switch = 3;
+        const int min_its_after_switch = this->param_.min_iter_after_switch_ms_wells_;
         // We also want to restrict the number of status switches to avoid oscillation between STOP<->OPEN
         const int max_status_switch = this->param_.max_well_status_switch_inner_iter_;
         int its_since_last_switch = min_its_after_switch;

opm/simulators/wells/StandardWellEval.cpp

@@ -157,6 +157,9 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
                    const std::vector<Scalar>& B_avg,
                    const Scalar maxResidualAllowed,
                    const Scalar tol_wells,
+                   const Scalar tol_bhp_eq,
+                   const Scalar tol_thp_eq,
+                   const Scalar tol_grup_eq,
                    const Scalar relaxed_tolerance_flow,
                    const bool relax_tolerance,
                    const bool well_is_stopped,
@@ -208,7 +211,7 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
 
     WellConvergence(baseif_).
         checkConvergenceControlEq(well_state,
-                                  {1.e3, 1.e4, 1.e-4, 1.e-6, maxResidualAllowed},
+                                  {tol_bhp_eq, tol_thp_eq, tol_wells, tol_grup_eq, maxResidualAllowed},
                                   std::abs(this->linSys_.residual()[0][Bhp]),
                                   well_is_stopped,
                                   report,

opm/simulators/wells/StandardWellEval.hpp

@@ -101,6 +101,9 @@ class StandardWellEval
                                          const std::vector<Scalar>& B_avg,
                                          const Scalar maxResidualAllowed,
                                          const Scalar tol_wells,
+                                         const Scalar tol_bhp_eq,
+                                         const Scalar tol_thp_eq,
+                                         const Scalar tol_grup_eq,
                                          const Scalar relaxed_tolerance_flow,
                                          const bool relax_tolerance,
                                          const bool well_is_stopped,

opm/simulators/wells/StandardWell_impl.hpp

@@ -1238,9 +1238,9 @@ namespace Opm
         auto& deferred_logger = groupStateHelper.deferredLogger();
         Scalar tol_wells = this->param_.tolerance_wells_;
         // use stricter tolerance for stopped wells and wells under zero rate target control.
-        constexpr Scalar stopped_factor = 1.e-4;
+        const Scalar stopped_factor = this->param_.tolerance_wells_stopped_factor_;
         // use stricter tolerance for dynamic thp to ameliorate network convergence
-        constexpr Scalar dynamic_thp_factor = 1.e-1;
+        const Scalar dynamic_thp_factor = this->param_.tolerance_wells_dynamic_thp_factor_;
         if (this->stoppedOrZeroRateTarget(groupStateHelper)) {
             tol_wells = tol_wells*stopped_factor;
         } else if (this->getDynamicThpLimit()) {
@@ -1252,6 +1252,9 @@ namespace Opm
                                                                          B_avg,
                                                                          this->param_.max_residual_allowed_,
                                                                          tol_wells,
+                                                                         this->param_.tolerance_well_bhp_eq_,
+                                                                         this->param_.tolerance_well_thp_eq_,
+                                                                         this->param_.tolerance_well_grup_eq_,
                                                                          this->param_.relaxed_tolerance_flow_well_,
                                                                          relax_tolerance,
                                                                          this->wellIsStopped(),
@@ -2479,7 +2482,7 @@ namespace Opm
         // Always take a few (more than one) iterations after a switch before allowing a new switch
         // The optimal number here is subject to further investigation, but it has been observerved
         // that unless this number is >1, we may get stuck in a cycle
-        constexpr int min_its_after_switch = 4;
+        const int min_its_after_switch = this->param_.min_iter_after_switch_wells_;
         // We also want to restrict the number of status switches to avoid oscillation between STOP<->OPEN
         const int max_status_switch = this->param_.max_well_status_switch_inner_iter_;
         int its_since_last_switch = min_its_after_switch;