Let the Newton relaxation factor recover after an oscillation (optional)

opm/simulators/flow/NonlinearSolver.cpp

@@ -146,6 +146,7 @@ NonlinearSolverParameters()
 
     // overload with given parameters
     relaxMax_ = Parameters::Get<Parameters::NewtonMaxRelax<Scalar>>();
+    relaxRecovery_ = Parameters::Get<Parameters::NewtonRelaxRecovery>();
 
     const auto& relaxationTypeString = Parameters::Get<Parameters::NewtonRelaxationType>();
     if (relaxationTypeString == "dampen") {
@@ -167,6 +168,7 @@ reset()
     relaxMax_ = 0.5;
     relaxIncrement_ = 0.1;
     relaxRelTol_ = 0.2;
+    relaxRecovery_ = 0;
 }
 
 template<class Scalar>
@@ -175,6 +177,11 @@ registerParameters()
 {
     Parameters::Register<Parameters::NewtonMaxRelax<Scalar>>
         ("The maximum relaxation factor of a Newton iteration");
+    Parameters::Register<Parameters::NewtonRelaxRecovery>
+        ("Number of consecutive Newton iterations with contracting residuals after which "
+         "the relaxation factor is walked back up towards 1. The factor otherwise only "
+         "ever decreases within a substep, so one oscillation early on damps every later "
+         "iteration. 0 (default) keeps the factor latched once reduced");
     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,9 @@ namespace Opm::Parameters {
 template<class Scalar>
 struct NewtonMaxRelax { static constexpr Scalar value = 0.5; };
 
+struct NewtonRelaxRecovery { static constexpr int value = 0; };
+
+
 struct NewtonRelaxationType { static constexpr auto value = "dampen"; };
 
 } // namespace Opm::Parameters
@@ -82,6 +85,7 @@ struct NonlinearSolverParameters
     Scalar relaxMax_;
     Scalar relaxIncrement_;
     Scalar relaxRelTol_;
+    int relaxRecovery_;
 
     NonlinearSolverParameters();
 
@@ -248,6 +252,11 @@ struct NonlinearSolverParameters
         }
 
         /// The greatest relaxation factor (i.e. smallest factor) allowed.
+        //! \brief Contracting iterations required before the relaxation factor recovers
+        //! (0 = never; the factor stays latched at its reduced value for the substep).
+        int relaxRecovery() const
+        { return param_.relaxRecovery_; }
+
         Scalar relaxMax() const
         { return param_.relaxMax_; }
 

opm/simulators/flow/NonlinearSystem.hpp

@@ -170,6 +170,9 @@ class NonlinearSystem
     ComponentName compNames_{};
     std::vector<std::vector<Scalar>> residual_norms_history_;
     Scalar current_relaxation_;
+    //! Consecutive iterations with contracting residuals, used to let
+    //! current_relaxation_ recover after an oscillation has passed.
+    int contracting_iterations_{0};
     GlobalEqVector dx_old_;
 };
 

opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp

@@ -210,6 +210,7 @@ nonlinearIteration(const SimulatorTimerInterface& timer,
         this->residual_norms_history_.clear();
         this->conv_monitor_.reset();
         this->current_relaxation_ = 1.0;
+        this->contracting_iterations_ = 0;
         this->dx_old_ = 0.0;
         this->convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
         this->convergence_reports_.back().report.reserve(11);
@@ -290,11 +291,40 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer,
             if (isOscillate) {
                 this->current_relaxation_ -= nonlinear_solver.relaxIncrement();
                 this->current_relaxation_ = std::max(this->current_relaxation_, nonlinear_solver.relaxMax());
+                this->contracting_iterations_ = 0;
                 if (this->terminalOutputEnabled()) {
                     OpmLog::info("    Oscillating behavior detected: Relaxation set to "
                                  + std::to_string(this->current_relaxation_));
                 }
             }
+            else if ((nonlinear_solver.relaxRecovery() > 0) &&
+                     (this->current_relaxation_ < 1.0))
+            {
+                // The factor is otherwise only ever reduced within a substep, so a single
+                // oscillation early on damps every later iteration. Let it recover once
+                // the residuals have contracted for a few iterations in a row.
+                const auto& hist = this->residual_norms_history_;
+                bool contracting = false;
+                if (hist.size() >= 2) {
+                    const auto normOf = [](const std::vector<Scalar>& norms)
+                    {
+                        return *std::max_element(norms.begin(), norms.end());
+                    };
+                    contracting = normOf(hist.back()) < normOf(hist[hist.size() - 2]);
+                }
+                this->contracting_iterations_ = contracting
+                    ? this->contracting_iterations_ + 1 : 0;
+                if (this->contracting_iterations_ >= nonlinear_solver.relaxRecovery()) {
+                    this->current_relaxation_ =
+                        std::min(Scalar{1.0},
+                                 this->current_relaxation_ + nonlinear_solver.relaxIncrement());
+                    this->contracting_iterations_ = 0;
+                    if (this->terminalOutputEnabled()) {
+                        OpmLog::info("    Residuals contracting: Relaxation restored to "
+                                     + std::to_string(this->current_relaxation_));
+                    }
+                }
+            }
             nonlinear_solver.stabilizeNonlinearUpdate(x, this->dx_old_, this->current_relaxation_);
         }