General controller

opm/simulators/flow/NonlinearSolver.hpp

@@ -36,6 +36,8 @@
 
 #include <opm/simulators/timestepping/SimulatorReport.hpp>
 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
+#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
+#include <opm/simulators/timestepping/TimeStepControl.hpp>
 
 #include <memory>
 
@@ -83,6 +85,7 @@ void registerNonlinearParameters();
     class NonlinearSolver
     {
         using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        using TimeStepper = AdaptiveTimeStepping<TypeTag>;
 
     public:
         // Solver parameters controlling nonlinear process.
@@ -161,7 +164,7 @@ void registerNonlinearParameters();
         }
 
 
-        SimulatorReportSingle step(const SimulatorTimerInterface& timer)
+        SimulatorReportSingle step(const SimulatorTimerInterface& timer, const TimeStepControlInterface& timeStepControl)
         {
             SimulatorReportSingle report;
             report.global_time = timer.simulationTimeElapsed();
@@ -208,6 +211,14 @@ void registerNonlinearParameters();
                 OPM_THROW_NOLOG(TooManyIterations, msg);
             }
 
+            if (!timeStepControl.timeStepAccepted(model_->relativeChange())) {
+                report.converged = false;
+                failureReport_ = report;
+
+                std::string msg = "Time step too large - Failed to satisfy the tolerance test.";
+                OPM_THROW_NOLOG(TimeSteppingBreakdown, msg);
+            }
+
             // Do model-specific post-step actions.
             report += model_->afterStep(timer);
             report.converged = true;

opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp

@@ -356,7 +356,7 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
             simulator_.problem().setSimulationReport(report_);
         } else {
             // solve for complete report step
-            auto stepReport = solver_->step(timer);
+            auto stepReport = solver_->step(timer, adaptiveTimeStepping_->timeStepControl());
             report_ += stepReport;
             if (terminalOutput_) {
                 std::ostringstream ss;
@@ -443,6 +443,9 @@ class SimulatorFullyImplicitBlackoil : private SerializableSim
     const Model& model() const
     { return solver_->model(); }
 
+    const TimeStepper& timeStepper() const
+    { return adaptiveTimeStepping_; }
+
 protected:
     //! \brief Load simulator state from hdf5 serializer.
     void loadState([[maybe_unused]] HDF5Serializer& serializer,

opm/simulators/timestepping/AdaptiveTimeStepping.cpp

@@ -136,6 +136,12 @@ void registerAdaptiveParameters()
     Parameters::Register<Parameters::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::Register<Parameters::TimeStepChoppingFactor>
+        ("Safety factor in the formula for the time step chopping after a time "
+        "step has failed to satisfy the tolerance criterion due to a too large time step");
+    Parameters::Register<Parameters::TimeStepControlSafetyFactor>
+        ("Safety factor to be multiplied with the tolerance parameter in the controllers "
+        "used to propose time steps");
 }
 
 } // namespace detail

opm/simulators/timestepping/AdaptiveTimeStepping.hpp

@@ -34,6 +34,8 @@
 
 #include <opm/simulators/utils/phaseUsageFromDeck.hpp>
 
+#include <boost/date_time/posix_time/posix_time.hpp>
+
 #include <fmt/format.h>
 
 #include <algorithm>
@@ -66,6 +68,8 @@ struct TimeStepControlGrowthDampingFactor { static constexpr double value = 3.2;
 struct TimeStepControlFileName { static constexpr auto value = "timesteps"; };
 struct MinTimeStepBeforeShuttingProblematicWellsInDays { static constexpr double value = 0.01; };
 struct MinTimeStepBasedOnNewtonIterations { static constexpr double value = 0.0; };
+struct TimeStepChoppingFactor { static constexpr double value = 0.8; };
+struct TimeStepControlSafetyFactor { static constexpr double value = 0.8; };
 
 } // namespace Opm::Parameters
 
@@ -232,8 +236,9 @@ void registerAdaptiveParameters();
 
                 SimulatorReportSingle substepReport;
                 std::string causeOfFailure;
+                bool tooLargeTimeStep = false;
                 try {
-                    substepReport = solver.step(substepTimer);
+                    substepReport = solver.step(substepTimer, *timeStepControl_);
 
                     if (solverVerbose_) {
                         // report number of linear iterations
@@ -258,6 +263,13 @@ void registerAdaptiveParameters();
                     logException_(e, solverVerbose_);
                     // since linearIterations is < 0 this will restart the solver
                 }
+                catch (const TimeSteppingBreakdown& e) {
+                    tooLargeTimeStep = true;
+                    substepReport = solver.failureReport();
+                    causeOfFailure = "Time step was too large";
+
+                    logException_(e, solverVerbose_);
+                }
                 catch (const NumericalProblem& e) {
                     substepReport = solver.failureReport();
                     causeOfFailure = "Solver convergence failure - Numerical problem encountered";
@@ -304,19 +316,26 @@ void registerAdaptiveParameters();
                     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) {
+
+                    double RC = relativeChange.relativeChange();
+                    if (solver.model().simulator().gridView().comm().rank() == 0)
+                    {
+                        std::cout << "RC: " << RC << std::endl;
+                        std::cout << "T: " << boost::posix_time::to_iso_extended_string(substepTimer.currentDateTime()) << std::endl;
+                    }
 
                     // 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;
                     double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
-                                                                               substepTimer.simulationTimeElapsed());
+                                                                               substepTimer);
 
                     assert(dtEstimate > 0);
                     // limit the growth of the timestep size by the growth factor
@@ -365,6 +384,57 @@ void registerAdaptiveParameters();
                     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 newTimeStep = timeStepChoppingFactor_ * 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);
 
@@ -481,6 +551,9 @@ void registerAdaptiveParameters();
         void setSuggestedNextStep(const double x)
         { suggestedNextTimestep_ = x; }
 
+        const TimeStepControlInterface& timeStepControl() const
+        { return *timeStepControl_; }
+
         void updateTUNING(double max_next_tstep, const Tuning& tuning)
         {
             restartFactor_ = tuning.TSFCNV;
@@ -516,6 +589,8 @@ void registerAdaptiveParameters();
             case TimeStepControlType::PID:
                 allocAndSerialize<PIDTimeStepControl>(serializer);
                 break;
+            case TimeStepControlType::General3rdOrder:
+                allocAndSerialize<General3rdOrderController>(serializer);
             }
             serializer(restartFactor_);
             serializer(growthFactor_);
@@ -553,6 +628,11 @@ void registerAdaptiveParameters();
             return serializationTestObject_<SimpleIterationCountTimeStepControl>();
         }
 
+        static AdaptiveTimeStepping<TypeTag> serializationTestObjectGeneral3rdOrder()
+        {
+            return serializationTestObject_<General3rdOrderController>();
+        }
+
         bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs)
         {
             if (timeStepControlType_ != rhs.timeStepControlType_ ||
@@ -575,6 +655,9 @@ void registerAdaptiveParameters();
             case TimeStepControlType::PID:
                 result = castAndComp<PIDTimeStepControl>(rhs);
                 break;
+            case TimeStepControlType::General3rdOrder:
+                result = castAndComp<General3rdOrderController>(rhs);
+                break;
             }
 
             return result &&
@@ -640,16 +723,19 @@ void registerAdaptiveParameters();
             // valid are "pid" and "pid+iteration"
             std::string control = Parameters::Get<Parameters::TimeStepControl>(); // "pid"
 
-            const double tol =  Parameters::Get<Parameters::TimeStepControlTolerance>(); // 1e-1
+            timeStepControlTolerance_ =  Parameters::Get<Parameters::TimeStepControlTolerance>(); // 1e-1
+            timeStepChoppingFactor_ = Parameters::Get<Parameters::TimeStepChoppingFactor>(); // 0.8
+            timeStepControlSafetyFactor_ = Parameters::Get<Parameters::TimeStepControlSafetyFactor>(); // 0.8
+
             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<Parameters::TimeStepControlTargetIterations>(); // 30
                 const double decayDampingFactor = Parameters::Get<Parameters::TimeStepControlDecayDampingFactor>(); // 1.0
                 const double growthDampingFactor = Parameters::Get<Parameters::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") {
@@ -660,7 +746,7 @@ void registerAdaptiveParameters();
                 // 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;
             }
@@ -684,6 +770,10 @@ void registerAdaptiveParameters();
                 timeStepControl_ = std::make_unique<HardcodedTimeStepControl>(filename);
                 timeStepControlType_ = TimeStepControlType::HardCodedTimeStep;
             }
+            else if (control == "general3rdorder") {
+                timeStepControl_ = std::make_unique<General3rdOrderController>(timeStepControlTolerance_, timeStepControlSafetyFactor_);
+                timeStepControlType_ = TimeStepControlType::General3rdOrder;
+            }
             else
                 OPM_THROW(std::runtime_error,
                           "Unsupported time step control selected " + control);
@@ -696,6 +786,7 @@ void registerAdaptiveParameters();
 
         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
@@ -710,6 +801,8 @@ void registerAdaptiveParameters();
         double timestepAfterEvent_;          //!< suggested size of timestep after an event
         bool useNewtonIteration_;            //!< use newton iteration count for adaptive time step control
         double minTimeStepBeforeShuttingProblematicWells_; //! < shut problematic wells when time step size in days are less than this
+        double timeStepChoppingFactor_;      //!< multiplicative factor for time step chopping after tolerance test fail
+        double timeStepControlSafetyFactor_; //!< factor multiplied with tolerance parameter in the adaptive time stepping controllers
     };
 }