Merge pull request #6280 from erikhide/general-controller-further-updates

Further updates and fixes to general controller

opm/simulators/flow/NonlinearSolver.hpp

@@ -176,8 +176,9 @@ struct NonlinearSolverParameters
                 OPM_THROW_NOLOG(TooManyIterations, msg);
             }
             auto relativeChange = model_->relativeChange();
-            if (timeStepControl && !timeStepControl->timeStepAccepted(relativeChange)) {
+            if (timeStepControl && !timeStepControl->timeStepAccepted(relativeChange, timer.currentStepLength())) {
                 report.converged = false;
+                report.time_step_rejected = true;
                 failureReport_ = report;
 
                 std::string msg = "Relative change in solution for time step was " + std::to_string(relativeChange) + ", which is larger than the tolerance accepted by the timestepping algorithm.";

opm/simulators/timestepping/AdaptiveTimeStepping.cpp

@@ -142,6 +142,19 @@ void registerAdaptiveParameters()
     Parameters::Register<Parameters::TimeStepControlRejectCompletedStep>
         ("(Only applicable for the general 3rd order controller.) Include rejection of completed "
          "time steps if the relative change is larger than the time step control tolerance");
+    Parameters::Register<Parameters::TimeStepControlToleranceTestVersion>
+        ("(Only applicable for the general 3rd order controller.) Ways to decide if the time step "
+         "should be rejected. Options: 'standard' and 'control-error-filtering'. The standard "
+         "version compares relative change to tolerance directly to decide if the time step should "
+         "be rejected, while the control-error-filtering version compares the relative change "
+         "in time step size to the max-reduction-time-step parameter.");
+    Parameters::Register<Parameters::TimeStepControlMaxReductionTimeStep>
+        ("(Only applicable for the general 3rd order controller, using 'control-error-filtering' "
+         "as time-step-control-tolerance-test-version) If the (proposed) relative change in time "
+         "step size is larger than this parameter, the time step will be rejected.");
+    Parameters::Register<Parameters::TimeStepControlParameters>
+        ("(Only applicable for the general 3rd order controller.) Parameters for the general "
+         "3rd order controller. Should be given as 'beta_1;beta_2;beta_3;alpha_2;alpha_3'.");
 }
 
 std::tuple<TimeStepControlType, std::unique_ptr<TimeStepControlInterface>, bool>
@@ -220,13 +233,21 @@ createController(const UnitSystem& unitSystem)
          }},
         {"general3rdorder",
          [tol]() {
-             const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>(); // 0.8
-             const bool rejectCompletedStep = Parameters::Get<Parameters::TimeStepControlRejectCompletedStep>(); // false
+             const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>();
+             const bool rejectCompletedStep = Parameters::Get<Parameters::TimeStepControlRejectCompletedStep>();
+             const std::string toleranceTestVersion = Parameters::Get<Parameters::TimeStepControlToleranceTestVersion>();
+             const double maxReductionTimeStep = Parameters::Get<Parameters::TimeStepControlMaxReductionTimeStep>();
+             const std::string parameters = Parameters::Get<Parameters::TimeStepControlParameters>();
+             const bool verbose = Parameters::Get<Parameters::TimeStepVerbosity>();
              return RetVal{
                  TimeStepControlType::General3rdOrder,
                  std::make_unique<General3rdOrderController>(tol,
                                                              safetyFactor,
-                                                             rejectCompletedStep),
+                                                             rejectCompletedStep,
+                                                             toleranceTestVersion,
+                                                             maxReductionTimeStep,
+                                                             parameters,
+                                                             verbose),
                  false
              };
         }},

opm/simulators/timestepping/AdaptiveTimeStepping.hpp

@@ -49,6 +49,9 @@ struct MinTimeStepBeforeShuttingProblematicWellsInDays { static constexpr double
 struct MinTimeStepBasedOnNewtonIterations { static constexpr double value = 0.0; };
 struct TimeStepControlSafetyFactor { static constexpr double value = 0.8; };
 struct TimeStepControlRejectCompletedStep { static constexpr bool value = false; };
+struct TimeStepControlToleranceTestVersion { static constexpr auto value = "standard"; };
+struct TimeStepControlMaxReductionTimeStep { static constexpr double value = 0.1; };
+struct TimeStepControlParameters { static constexpr auto value = "0.125;0.25;0.125;0.75;0.25"; };
 
 } // namespace Opm::Parameters
 

opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp

@@ -825,7 +825,15 @@ run()
             this->substep_timer_.setLastStepFailed(true);
             checkTimeStepMaxRestartLimit_(restarts);
 
-            const double new_time_step = restartFactor_() * dt;
+            double new_time_step = restartFactor_() * dt;
+            if (substep_report.time_step_rejected) {
+                const double tol = Parameters::Get<Parameters::TimeStepControlTolerance>();
+                const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>();
+                const double temp_time_step = std::sqrt(safetyFactor * tol / solver_().model().relativeChange()) * dt;
+                if (temp_time_step < dt) { // added in case suggested time step is not a reduction
+                    new_time_step = temp_time_step;
+                }
+            }
             checkTimeStepMinLimit_(new_time_step);
             bool wells_shut = false;
             if (new_time_step > minTimeStepBeforeClosingWells_()) {

opm/simulators/timestepping/SimulatorReport.cpp

@@ -35,7 +35,7 @@ namespace Opm
         return SimulatorReportSingle{1.0, 2.0, 3.0, 4.0, 5.0, 6.0,
                                      7.0, 8.0, 9.0, 10.0, 11.0, 12.0,
                                      13, 14, 15, 16, 17, 18,
-                                     true, false, 19, 20.0, 21.0,
+                                     true, false, false, 19, 20.0, 21.0,
                                      22, 23, 24, 25, 26, 27, 28};
     }
 
@@ -60,6 +60,7 @@ namespace Opm
                this->min_linear_iterations == rhs.min_linear_iterations &&
                this->max_linear_iterations == rhs.max_linear_iterations &&
                this->converged == rhs.converged &&
+               this->time_step_rejected == rhs.time_step_rejected &&
                this->well_group_control_changed == rhs.well_group_control_changed &&
                this->exit_status == rhs.exit_status &&
                this->global_time == rhs.global_time &&

opm/simulators/timestepping/SimulatorReport.hpp

@@ -53,6 +53,7 @@ namespace Opm
         unsigned int max_linear_iterations = 0;
 
         bool converged = false;
+        bool time_step_rejected = false;
         bool well_group_control_changed = false;
         int exit_status = EXIT_SUCCESS;
 
@@ -101,6 +102,7 @@ namespace Opm
             serializer(min_linear_iterations);
             serializer(max_linear_iterations);
             serializer(converged);
+            serializer(time_step_rejected);
             serializer(well_group_control_changed);
             serializer(exit_status);
             serializer(global_time);

opm/simulators/timestepping/TimeStepControl.cpp

@@ -28,6 +28,7 @@
 #include <string>
 #include <fstream>
 #include <iostream>
+#include <sstream>
 #include <limits>
 
 #include <opm/common/ErrorMacros.hpp>
@@ -292,83 +293,150 @@ namespace Opm
     General3rdOrderController::General3rdOrderController( const double tolerance,
                                                           const double safetyFactor,
                                                           const bool rejectCompletedStep,
+                                                          const std::string& toleranceTestVersion,
+                                                          const double maxReductionTimeStep,
+                                                          const std::string& parameters,
                                                           const bool verbose)
         : tolerance_( tolerance )
         , safetyFactor_( safetyFactor )
         , rejectCompletedStep_( rejectCompletedStep )
-        , errors_( 3, tolerance_ )
-        , timeSteps_ ( 3, 1.0 )
+        , maxReductionTimeStep_( maxReductionTimeStep )
         , verbose_( verbose )
-    {}
+    {
+        std::istringstream parameters_stream(parameters);
+        std::string value;
+        int counter = 0;
+        try {
+            while (std::getline(parameters_stream, value, ';'))
+            {
+                if (counter < 3) {
+                    beta_[counter] = std::stod(value);
+                }
+                else {
+                    alpha_[counter - 3] = std::stod(value);
+                }
+                counter++;
+            }
+        }
+        catch (const std::exception& e) {
+            OPM_THROW(std::runtime_error, "The parameter string for the time step controller is invalid: " + parameters);
+        }
+
+        if (toleranceTestVersion == "standard") {
+            toleranceTestVersion_ = ToleranceTestVersions::Standard;
+        }
+        else if (toleranceTestVersion == "control-error-filtering") {
+            toleranceTestVersion_ = ToleranceTestVersions::ControlErrorFiltering;
+        }
+        else {
+            OPM_THROW(std::runtime_error, "Unsupported tolerance test version given: " + toleranceTestVersion);
+        }
+    }
 
     General3rdOrderController
     General3rdOrderController::serializationTestObject()
     {
-        General3rdOrderController result(1.0, 0.8, true);
+        General3rdOrderController result(1.0, 2.0, false, "standard", 3.0, "0.125;0.25;0.125;0.75;0.25", false);
         result.errors_ = {2.0, 3.0};
 
         return result;
     }
 
     double General3rdOrderController::
-    computeTimeStepSize(const double dt, const int /*iterations */, const RelativeChangeInterface& relChange, const AdaptiveSimulatorTimer& substepTimer) const
+    computeTimeStepSize(const double dt, const int /* iterations */, const RelativeChangeInterface& /* relChange */, const AdaptiveSimulatorTimer& substepTimer) const
     {
-        // Shift errors and time steps
-        for( int i = 0; i < 2; ++i )
+        if (errors_[0] == 0 || errors_[1] == 0 || errors_[2] == 0)
         {
-            errors_[i] = errors_[i+1];
-            timeSteps_[i] = timeSteps_[i+1];
+            if ( verbose_ )
+            {
+                OpmLog::info("The solution between time steps does not change, there is no time step constraint from the controller.");
+            }
+            return std::numeric_limits<double>::max();
         }
 
-        // Store new error and time step
-        const double error = relChange.relativeChange();
-        errors_[2] = error;
-        timeSteps_[2] = dt;
         for( int i = 0; i < 2; ++i )
         {
             assert(std::isfinite(errors_[i]));
         }
 
-        if (errors_[0] == 0 || errors_[1] == 0 || errors_[2] == 0.)
+        // Use an I controller after report time steps
+        if (substepTimer.currentStepNum() < 3)
         {
-            if ( verbose_ )
-                OpmLog::info("The solution between time steps does not change, there is no time step constraint from the controller.");
-            return std::numeric_limits<double>::max();
-        }
-        // Use an I controller after report time steps or chopped time steps
-        else if (substepTimer.currentStepNum() < 3 || substepTimer.lastStepFailed() || counterSinceFailure_ > 0)
-        {
-            if (substepTimer.lastStepFailed() || counterSinceFailure_ > 0)
-                counterSinceFailure_++;
-            if (counterSinceFailure_ > 1)
-                counterSinceFailure_ = 0;
-            const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], 0.35);
-            if( verbose_ )
+            controllerVersion_ = InternalControlVersions::IController;
+            const double newDt = dt * timeStepFactor(errors_, timeSteps_);
+            if (verbose_)
+            {
                 OpmLog::info(fmt::format("Computed step size (pow): {} days", unit::convert::to( newDt, unit::day )));
+            }
             return newDt;
         }
         // Use the general third order controller for all other time steps
         else
         {
-            const std::array<double, 3> beta = { 0.125, 0.25, 0.125 };
-            const std::array<double, 2> alpha = { 0.375, 0.125 };
-            const double newDt = dt * std::pow(safetyFactor_ * tolerance_ / errors_[2], beta[0]) *
-                                      std::pow(safetyFactor_ * tolerance_ / errors_[1], beta[1]) *
-                                      std::pow(safetyFactor_ * tolerance_ / errors_[0], beta[2]) *
-                                      std::pow(timeSteps_[2] / timeSteps_[1], -alpha[0]) *
-                                      std::pow(timeSteps_[1] / timeSteps_[0], -alpha[1]);
-            if( verbose_ )
+            controllerVersion_ = InternalControlVersions::General3rdOrder;
+            const double newDt = dt * timeStepFactor(errors_, timeSteps_);
+            if (verbose_)
+            {
                 OpmLog::info(fmt::format("Computed step size (pow): {} days", unit::convert::to( newDt, unit::day )));
+            }
             return newDt;
         }
     }
 
-    bool General3rdOrderController::
-    timeStepAccepted(const double error) const
+    double General3rdOrderController::
+    timeStepFactor(const std::array<double, 3>& errors, const std::array<double, 3>& timeSteps) const
     {
-        if (rejectCompletedStep_ && error > tolerance_)
-            return false;
-        return true;
+        if (controllerVersion_ == InternalControlVersions::General3rdOrder)
+        {
+            return std::pow(safetyFactor_ * tolerance_ / errors[2], beta_[0]) *
+                   std::pow(safetyFactor_ * tolerance_ / errors[1], beta_[1]) *
+                   std::pow(safetyFactor_ * tolerance_ / errors[0], beta_[2]) *
+                   std::pow(timeSteps[2] / timeSteps[1], -alpha_[0]) *
+                   std::pow(timeSteps[1] / timeSteps[0], -alpha_[1]);
+        }
+        // controllerVersion_ == InternalControlVersions::IController
+        return std::pow(safetyFactor_ * tolerance_ / errors[2], 0.35);
+    }
+
+    bool General3rdOrderController::
+    timeStepAccepted(const double error, const double timeStepJustCompleted) const
+    {
+        bool acceptTimeStep = true;
+
+        // Reject time step if chosen tolerance test version fails
+        if (toleranceTestVersion_ == ToleranceTestVersions::Standard)
+        {
+            if (rejectCompletedStep_ && error > tolerance_)
+            {
+                acceptTimeStep = false;
+            }
+        }
+        else if (toleranceTestVersion_ == ToleranceTestVersions::ControlErrorFiltering)
+        {
+            const std::array<double, 3> tempErrors{errors_[0], errors_[1], error};
+            const std::array<double, 3> tempTimeSteps{timeSteps_[0], timeSteps_[1], timeStepJustCompleted};
+            const double stepFactor = timeStepFactor(tempErrors, tempTimeSteps);
+            if (rejectCompletedStep_ && stepFactor < maxReductionTimeStep_)
+            {
+                acceptTimeStep = false;
+            }
+        }
+
+        if (acceptTimeStep)
+        {
+            // Shift errors and time steps
+            for( int i = 0; i < 2; ++i )
+            {
+                errors_[i] = errors_[i+1];
+                timeSteps_[i] = timeSteps_[i+1];
+            }
+
+            // Store new error and time step
+            errors_[2] = error;
+            timeSteps_[2] = timeStepJustCompleted;
+        }
+        
+        return acceptTimeStep;
     }
 
     bool General3rdOrderController::operator==(const General3rdOrderController& ctrl) const
@@ -378,6 +446,11 @@ namespace Opm
                this->rejectCompletedStep_ == ctrl.rejectCompletedStep_ &&
                this->errors_ == ctrl.errors_ &&
                this->timeSteps_ == ctrl.timeSteps_ &&
+               this->beta_ == ctrl.beta_ &&
+               this->alpha_ == ctrl.alpha_ &&
+               this->controllerVersion_ == ctrl.controllerVersion_ &&
+               this->toleranceTestVersion_ == ctrl.toleranceTestVersion_ &&
+               this->maxReductionTimeStep_ == ctrl.maxReductionTimeStep_ &&
                this->verbose_ == ctrl.verbose_;
     }
 

opm/simulators/timestepping/TimeStepControl.hpp

@@ -37,6 +37,16 @@ namespace Opm
       General3rdOrder,
     };
 
+    enum class ToleranceTestVersions {
+        Standard,
+        ControlErrorFiltering,
+    };
+
+    enum class InternalControlVersions {
+        IController,
+        General3rdOrder,
+    };
+
     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ///
     ///  A simple iteration count based adaptive time step control.
@@ -66,7 +76,7 @@ namespace Opm
                                    const RelativeChangeInterface& /* relativeChange */,
                                    const AdaptiveSimulatorTimer& /* substepTimer */ ) const override;
         
-        bool timeStepAccepted(const double /* error */) const override { return true; }
+        bool timeStepAccepted(const double /* error */, const double /* timeStepJustCompleted */) const override { return true; }
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -119,7 +129,7 @@ namespace Opm
                                    const RelativeChangeInterface& relativeChange,
                                    const AdaptiveSimulatorTimer& /* substepTimer */ ) const override;
 
-        bool timeStepAccepted(const double /* error */) const override { return true; }
+        bool timeStepAccepted(const double /* error */, const double /* timeStepJustCompleted */) const override { return true; }
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -170,7 +180,7 @@ namespace Opm
                                    const RelativeChangeInterface& relativeChange,
                                    const AdaptiveSimulatorTimer& /* substepTimer */ ) const override;
 
-        bool timeStepAccepted(const double /* error */) const override { return true; }
+        bool timeStepAccepted(const double /* error */, const double /* timeStepJustCompleted */) const override { return true; }
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -201,19 +211,26 @@ namespace Opm
     public:
         static constexpr TimeStepControlType Type = TimeStepControlType::General3rdOrder;
 
-        General3rdOrderController( const double tolerance = 1e-3,
-                                   const double safetyFactor = 0.8,
-                                   const bool rejectCompletedStep = false,
-                                   const bool verbose = false );
+        General3rdOrderController() = default;
+
+        General3rdOrderController( const double tolerance,
+                                   const double safetyFactor,
+                                   const bool rejectCompletedStep,
+                                   const std::string& toleranceTestVersion,
+                                   const double maxReductionTimeStep,
+                                   const std::string& parameters,
+                                   const bool verbose);
 
         static General3rdOrderController serializationTestObject();
 
         double computeTimeStepSize(const double dt,
-                                   const int /*iterations */,
-                                   const RelativeChangeInterface& relativeChange,
+                                   const int /* iterations */,
+                                   const RelativeChangeInterface& /* relativeChange */,
                                    const AdaptiveSimulatorTimer& substepTimer) const override;
 
-        bool timeStepAccepted(const double error) const override;
+        double timeStepFactor(const std::array<double, 3>& errors, const std::array<double, 3>& timeSteps) const;
+
+        bool timeStepAccepted(const double error, const double timeStepJustCompleted) const override;
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -223,6 +240,11 @@ namespace Opm
             serializer(rejectCompletedStep_);
             serializer(errors_);
             serializer(timeSteps_);
+            serializer(beta_);
+            serializer(alpha_);
+            serializer(controllerVersion_);
+            serializer(toleranceTestVersion_);
+            serializer(maxReductionTimeStep_);
             serializer(verbose_);
         }
 
@@ -230,12 +252,16 @@ namespace Opm
 
 
     protected:
-        const double tolerance_ = 1e-3;
+        const double tolerance_ = 0.1;
         const double safetyFactor_ = 0.8;
         const bool rejectCompletedStep_ = false;
-        mutable std::vector<double> errors_{};
-        mutable std::vector<double> timeSteps_{};
-        mutable int counterSinceFailure_ = 0;
+        mutable std::array<double, 3> errors_{};
+        mutable std::array<double, 3> timeSteps_{};
+        mutable std::array<double, 3> beta_{0.125, 0.25, 0.125};
+        mutable std::array<double, 2> alpha_{0.75, 0.25};
+        mutable InternalControlVersions controllerVersion_{InternalControlVersions::IController};
+        ToleranceTestVersions toleranceTestVersion_{ToleranceTestVersions::Standard};
+        const double maxReductionTimeStep_ = 0.1;
         const bool verbose_ = false;
     };
 
@@ -266,7 +292,7 @@ namespace Opm
                                    const RelativeChangeInterface& /*relativeChange */,
                                    const AdaptiveSimulatorTimer& substepTimer) const override;
 
-        bool timeStepAccepted(const double /* error */) const override { return true; }
+        bool timeStepAccepted(const double /* error */, const double /* timeStepJustCompleted */) const override { return true; }
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)

opm/simulators/timestepping/TimeStepControlInterface.hpp

@@ -60,7 +60,9 @@ namespace Opm
         /// \return suggested time step size for the next step
         virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange , const AdaptiveSimulatorTimer& substepTimer) const = 0;
 
-        virtual bool timeStepAccepted(const double error) const = 0;
+        /// For the general 3rd order controller, the internal shifting of errors and time steps happens here, and
+        /// hence this method needs to be called for (after) each time step
+        virtual bool timeStepAccepted(const double error, const double timeStepJustCompleted) const = 0;
 
         /// virtual destructor (empty)
         virtual ~TimeStepControlInterface () {}