Include rejection of time step if tolerance test fails

opm/simulators/flow/NonlinearSolver.hpp

@@ -24,6 +24,7 @@
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
+#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
 
 #include <opm/models/utils/parametersystem.hh>
 #include <opm/models/utils/propertysystem.hh>
@@ -107,6 +108,7 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
     class NonlinearSolver
     {
         using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        using TimeStepper = AdaptiveTimeStepping<TypeTag>;
 
     public:
         // Solver parameters controlling nonlinear process.
@@ -142,6 +144,8 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
 
             static void registerParameters()
             {
+                TimeStepper::registerParameters();
+
                 Parameters::registerParam<TypeTag, Properties::NewtonMaxRelax>
                     ("The maximum relaxation factor of a Newton iteration");
                 Parameters::registerParam<TypeTag, Properties::NewtonMaxIterations>
@@ -239,6 +243,17 @@ void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld,
                 OPM_THROW_NOLOG(TooManyIterations, msg);
             }
 
+            double error = model_->relativeChange();
+            const double timeStepControlTolerance = Parameters::get<TypeTag, Properties::TimeStepControlTolerance>();
+
+            if (error > timeStepControlTolerance) {
+                report.converged = false;
+                failureReport_ = report;
+
+                std::string msg = "Time step too large - Failed to satisfy the tolerance test, since the error (" + std::to_string(error) + ") was larger than the tolerance (" + std::to_string(timeStepControlTolerance) + ").";
+                OPM_THROW_NOLOG(TimeSteppingBreakdown, msg);
+            }
+
             // Do model-specific post-step actions.
             report += model_->afterStep(timer);
             report.converged = true;

opm/simulators/timestepping/AdaptiveTimeStepping.hpp

@@ -123,6 +123,10 @@ template<class TypeTag, class MyTypeTag>
 struct MinTimeStepBasedOnNewtonIterations {
     using type = UndefinedProperty;
 };
+template<class TypeTag, class MyTypeTag>
+struct TimeStepSafetyFactor {
+    using type = UndefinedProperty;
+};
 
 template<class TypeTag>
 struct SolverContinueOnConvergenceFailure<TypeTag, TTag::FlowTimeSteppingParameters> {
@@ -201,6 +205,10 @@ struct MinTimeStepBasedOnNewtonIterations<TypeTag, TTag::FlowTimeSteppingParamet
     using type = GetPropType<TypeTag, Scalar>;
     static constexpr type value = 0.0;
 };
+template<class TypeTag>
+struct TimeStepSafetyFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
+    static constexpr double value = 0.8;
+};
 
 } // namespace Opm::Properties
 
@@ -354,6 +362,9 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
             Parameters::registerParam<TypeTag, Properties::MinTimeStepBasedOnNewtonIterations>
                 ("The minimum time step size (in days for field and metric unit and hours for lab unit) "
                  "can be reduced to based on newton iteration counts");
+            Parameters::registerParam<TypeTag, Properties::TimeStepSafetyFactor>
+                ("Safety factor in the formula for the time step cutting after a "
+                "time step has failed to satisfy the tolerance criterion");
         }
 
         /** \brief  step method that acts like the solver::step method
@@ -413,6 +424,7 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
 
                 SimulatorReportSingle substepReport;
                 std::string causeOfFailure;
+                bool tooLargeTimeStep = false;
                 try {
                     substepReport = solver.step(substepTimer);
 
@@ -435,6 +447,13 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
                     logException_(e, solverVerbose_);
                     // since linearIterations is < 0 this will restart the solver
                 }
+                catch (const TimeSteppingBreakdown& e) {
+                    tooLargeTimeStep = true;
+                    substepReport = solver.failureReport();
+                    causeOfFailure = "Error in time stepping exceeded tolerance";
+
+                    logException_(e, solverVerbose_);
+                }
                 catch (const NumericalProblem& e) {
                     substepReport = solver.failureReport();
                     causeOfFailure = "Solver convergence failure - Numerical problem encountered";
@@ -481,14 +500,14 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
                     OpmLog::problem(msg);
                 }
 
+                // create object to compute the time error, simply forwards the call to the model
+                SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
+
                 if (substepReport.converged || continue_on_uncoverged_solution) {
 
                     // advance by current dt
                     ++substepTimer;
 
-                    // create object to compute the time error, simply forwards the call to the model
-                    SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
-
                     // compute new time step estimate
                     const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
                         : substepReport.total_linear_iterations;
@@ -542,6 +561,58 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
                     substepTimer.setLastStepFailed(false);
 
                 }
+                else if (tooLargeTimeStep) {
+
+                    substepTimer.setLastStepFailed(true);
+
+                    // If we have restarted (i.e. cut the timestep) too
+                    // many times, we have failed and throw an exception.
+                    if (restarts >= solverRestartMax_) {
+                        const auto msg = fmt::format(
+                            "Solver failed to satisfy adaptive time stepping requirement."
+                        );
+                        if (solverVerbose_) {
+                            OpmLog::error(msg);
+                        }
+                        // Use throw directly to prevent file and line
+                        throw TimeSteppingBreakdown{msg};
+                    }
+
+                    // The new, chopped timestep.
+                    const double timeStepSafetyFactor = Parameters::get<TypeTag, Properties::TimeStepSafetyFactor>();
+                    const double newTimeStep = timeStepSafetyFactor * dt * timeStepControlTolerance_ / relativeChange.relativeChange();
+
+                    // If we have restarted (i.e. cut the timestep) too
+                    // much, we have failed and throw an exception.
+                    if (newTimeStep < minTimeStep_) {
+                        const auto msg = fmt::format(
+                            "Solver failed to converge after cutting timestep to {}\n"
+                            "which is the minimum threshold given by option --solver-min-time-step\n",
+                            minTimeStep_
+                        );
+                        if (solverVerbose_) {
+                            OpmLog::error(msg);
+                        }
+                        // Use throw directly to prevent file and line
+                        throw TimeSteppingBreakdown{msg};
+                    }
+
+                    // Define utility function for chopping timestep.
+                    auto chopTimestep = [&]() {
+                        substepTimer.provideTimeStepEstimate(newTimeStep);
+                        if (solverVerbose_) {
+                            const auto msg = fmt::format(
+                                "{}\nTimestep chopped to {} days\n",
+                                causeOfFailure,
+                                std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day))
+                            );
+                            OpmLog::problem(msg);
+                        }
+                        ++restarts;
+                    };
+
+                    chopTimestep();
+                }
                 else { // in case of no convergence
                     substepTimer.setLastStepFailed(true);
 
@@ -817,16 +888,16 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
             // valid are "pid" and "pid+iteration"
             std::string control = Parameters::get<TypeTag, Properties::TimeStepControl>(); // "pid"
 
-            const double tol =  Parameters::get<TypeTag, Properties::TimeStepControlTolerance>(); // 1e-1
+            timeStepControlTolerance_ =  Parameters::get<TypeTag, Properties::TimeStepControlTolerance>(); // 1e-1
             if (control == "pid") {
-                timeStepControl_ = std::make_unique<PIDTimeStepControl>(tol);
+                timeStepControl_ = std::make_unique<PIDTimeStepControl>(timeStepControlTolerance_);
                 timeStepControlType_ = TimeStepControlType::PID;
             }
             else if (control == "pid+iteration") {
                 const int iterations =  Parameters::get<TypeTag, Properties::TimeStepControlTargetIterations>(); // 30
                 const double decayDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlDecayDampingFactor>(); // 1.0
                 const double growthDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlGrowthDampingFactor>(); // 3.2
-                timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, tol);
+                timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, timeStepControlTolerance_);
                 timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
             }
             else if (control == "pid+newtoniteration") {
@@ -837,7 +908,7 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
                 // the min time step can be reduced by the newton iteration numbers
                 double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
                 timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor,
-                                                                                         growthDampingFactor, tol, minTimeStepReducedByIterations);
+                                                                                         growthDampingFactor, timeStepControlTolerance_, minTimeStepReducedByIterations);
                 timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
                 useNewtonIteration_ = true;
             }
@@ -873,6 +944,7 @@ std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr
 
         TimeStepControlType timeStepControlType_; //!< type of time step control object
         TimeStepController timeStepControl_; //!< time step control object
+        double timeStepControlTolerance_;    //!< tolerance for the adaptive time stepping
         double restartFactor_;               //!< factor to multiply time step with when solver fails to converge
         double growthFactor_;                //!< factor to multiply time step when solver recovered from failed convergence
         double maxGrowth_;                   //!< factor that limits the maximum growth of a time step

opm/simulators/timestepping/TimeStepControl.cpp

@@ -184,15 +184,7 @@ namespace Opm
             assert(std::isfinite(errors_[i]));
         }
 
-        if( errors_[2] > tol_ )
-        {
-            // adjust dt by given tolerance
-            const double newDt = dt * tol_ / error;
-            if ( verbose_ )
-                    OpmLog::info(fmt::format("Computed step size (tol): {} days", unit::convert::to( newDt, unit::day )));
-            return newDt;
-        }
-        else if (errors_[1] == 0 || errors_[2] == 0.)
+        if (errors_[0] == 0 || errors_[1] == 0 || errors_[2] == 0.)
         {
             if ( verbose_ )
                 OpmLog::info("The solution between time steps does not change, there is no time step constraint from the PID time step control ");
@@ -206,7 +198,7 @@ namespace Opm
             const double kD = 0.01 ;
             const double newDt = (dt * std::pow( errors_[ 1 ] / errors_[ 2 ], kP ) *
                                  std::pow( tol_         / errors_[ 2 ], kI ) *
-                                 std::pow( errors_[0]*errors_[0]/errors_[ 1 ]/errors_[ 2 ], kD ));
+                                 std::pow( errors_[1]*errors_[1]/errors_[ 0 ]/errors_[ 2 ], kD ));
             if( verbose_ )
                 OpmLog::info(fmt::format("Computed step size (pow): {} days", unit::convert::to( newDt, unit::day )));
             return newDt;