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

Further updates and fixes to general controller

Author: atgeirr

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;
         }
     }
 
+    double General3rdOrderController::
+    timeStepFactor(const std::array<double, 3>& errors, const std::array<double, 3>& timeSteps) const
+    {
+        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
+    timeStepAccepted(const double error, const double timeStepJustCompleted) const
     {
-        if (rejectCompletedStep_ && error > tolerance_)
-            return false;
-        return true;
+        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_;
     }