Merge pull request #5974 from erikhide/general-controller-updated

General third order time step controller

opm/simulators/timestepping/AdaptiveTimeStepping.cpp

@@ -136,6 +136,9 @@ 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::TimeStepControlSafetyFactor>
+        ("Value to be multiplied with the time step control tolerance to ensure that the target "
+         "relative change is lower than the tolerance");
 }
 
 std::tuple<TimeStepControlType, std::unique_ptr<TimeStepControlInterface>, bool>
@@ -212,6 +215,16 @@ createController(const UnitSystem& unitSystem)
                  true
              };
          }},
+        {"general3rdorder",
+         [tol]() {
+             const double safetyFactor = Parameters::Get<Parameters::TimeStepControlSafetyFactor>(); // 0.8
+             return RetVal{
+                 TimeStepControlType::General3rdOrder,
+                 std::make_unique<General3rdOrderController>(tol,
+                                                             safetyFactor),
+                 false
+             };
+        }},
         {"hardcoded",
          []() {
              const std::string filename = Parameters::Get<Parameters::TimeStepControlFileName>(); // "timesteps"

opm/simulators/timestepping/AdaptiveTimeStepping.hpp

@@ -57,6 +57,7 @@ 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 TimeStepControlSafetyFactor { static constexpr double value = 0.8; };
 
 } // namespace Opm::Parameters
 
@@ -173,7 +174,7 @@ private:
         const SimulatorTimer& simulatorTimer_() const;
         boost::posix_time::ptime startDateTime_() const;
         double timeStepControlComputeEstimate_(
-            const double dt, const int iterations, double elapsed) const;
+            const double dt, const int iterations, AdaptiveSimulatorTimer& substepTimer) const;
         bool timeStepVerbose_() const;
         void updateSuggestedNextStep_();
         bool useNewtonIteration_() const;
@@ -229,6 +230,7 @@ public:
     static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID();
     static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt();
     static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple();
+    static AdaptiveTimeStepping<TypeTag> serializationTestObject3rdOrder();
 
 private:
     void maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_step,

opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp

@@ -24,6 +24,7 @@
 #ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP
 #include <config.h>
 #include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
+#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
 #endif
 
 #include <opm/common/Exceptions.hpp>
@@ -144,6 +145,9 @@ operator==(const AdaptiveTimeStepping<TypeTag>& rhs)
     case TimeStepControlType::PID:
         result = castAndComp<PIDTimeStepControl>(rhs);
         break;
+    case TimeStepControlType::General3rdOrder:
+        result = castAndComp<General3rdOrderController>(rhs);
+        break;
     }
 
     return result &&
@@ -233,6 +237,9 @@ serializeOp(Serializer& serializer)
     case TimeStepControlType::PID:
         allocAndSerialize<PIDTimeStepControl>(serializer);
         break;
+    case TimeStepControlType::General3rdOrder:
+        allocAndSerialize<General3rdOrderController>(serializer);
+        break;
     }
     serializer(this->restart_factor_);
     serializer(this->growth_factor_);
@@ -282,6 +289,14 @@ serializationTestObjectSimple()
     return serializationTestObject_<SimpleIterationCountTimeStepControl>();
 }
 
+template<class TypeTag>
+AdaptiveTimeStepping<TypeTag>
+AdaptiveTimeStepping<TypeTag>::
+serializationTestObject3rdOrder()
+{
+    return serializationTestObject_<General3rdOrderController>();
+}
+
 
 template<class TypeTag>
 void
@@ -771,7 +786,7 @@ run()
 
             const int iterations = getNumIterations_(substep_report);
             auto dt_estimate = timeStepControlComputeEstimate_(
-                                     dt, iterations, this->substep_timer_.simulationTimeElapsed());
+                                     dt, iterations, this->substep_timer_);
 
             assert(dt_estimate > 0);
             dt_estimate = maybeRestrictTimeStepGrowth_(dt, dt_estimate, restarts);
@@ -1223,12 +1238,12 @@ template <class TypeTag>
 template <class Solver>
 double
 AdaptiveTimeStepping<TypeTag>::SubStepIteration<Solver>::
-timeStepControlComputeEstimate_(const double dt, const int iterations, double elapsed) const
+timeStepControlComputeEstimate_(const double dt, const int iterations, AdaptiveSimulatorTimer& substepTimer) const
 {
     // create object to compute the time error, simply forwards the call to the model
     SolutionTimeErrorSolverWrapper<Solver> relative_change{solver_()};
     return this->adaptive_time_stepping_.time_step_control_->computeTimeStepSize(
-        dt, iterations, relative_change, elapsed);
+        dt, iterations, relative_change, substepTimer);
 }
 
 template <class TypeTag>

opm/simulators/timestepping/TimeStepControl.cpp

@@ -20,6 +20,7 @@
 */
 #include <config.h>
 #include <algorithm>
+#include <array>
 #include <cassert>
 #include <cmath>
 #include <iostream>
@@ -33,6 +34,7 @@
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/input/eclipse/Units/Units.hpp>
 #include <opm/simulators/timestepping/TimeStepControl.hpp>
+#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
 
 #include <fmt/format.h>
 
@@ -73,7 +75,7 @@ namespace Opm
     }
 
     double SimpleIterationCountTimeStepControl::
-    computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const double /*simulationTimeElapsed */) const
+    computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& /* relativeChange */, const AdaptiveSimulatorTimer& /* substepTimer */) const
     {
         double dtEstimate = dt ;
 
@@ -134,10 +136,10 @@ namespace Opm
     }
 
     double HardcodedTimeStepControl::
-    computeTimeStepSize( const double /*dt */, const int /*iterations */, const RelativeChangeInterface& /* relativeChange */ , const double simulationTimeElapsed) const
+    computeTimeStepSize( const double /*dt */, const int /*iterations */, const RelativeChangeInterface& /* relativeChange */ , const AdaptiveSimulatorTimer& substepTimer) const
     {
-        auto nextTime = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), simulationTimeElapsed);
+        auto nextTime = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), substepTimer.simulationTimeElapsed());
-        return (*nextTime - simulationTimeElapsed);
+        return (*nextTime - substepTimer.simulationTimeElapsed());
     }
 
     bool HardcodedTimeStepControl::operator==(const HardcodedTimeStepControl& ctrl) const
@@ -170,7 +172,7 @@ namespace Opm
     }
 
     double PIDTimeStepControl::
-    computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relChange, const double /*simulationTimeElapsed */) const
+    computeTimeStepSize( const double dt, const int /* iterations */, const RelativeChangeInterface& relChange, const AdaptiveSimulatorTimer& /* substepTimer */) const
     {
         // shift errors
         for( int i=0; i<2; ++i ) {
@@ -249,9 +251,9 @@ namespace Opm
     }
 
     double PIDAndIterationCountTimeStepControl::
-    computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange,  const double simulationTimeElapsed ) const
+    computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange,  const AdaptiveSimulatorTimer& substepTimer) const
     {
-        double dtEstimatePID = PIDTimeStepControl :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed);
+        double dtEstimatePID = PIDTimeStepControl :: computeTimeStepSize( dt, iterations, relChange, substepTimer);
 
         // adjust timesteps based on target iteration
         double dtEstimateIter;
@@ -279,4 +281,93 @@ namespace Opm
                this->minTimeStepBasedOnIterations_ == ctrl.minTimeStepBasedOnIterations_;
     }
 
+
+
+    ////////////////////////////////////////////////////////////
+    //
+    //  General3rdOrderController  Implementation
+    //
+    ////////////////////////////////////////////////////////////
+
+    General3rdOrderController::General3rdOrderController( const double tolerance,
+                                                          const double safetyFactor,
+                                                          const bool verbose)
+        : tolerance_( tolerance )
+        , safetyFactor_( safetyFactor )
+        , errors_( 3, tolerance_ )
+        , timeSteps_ ( 3, 1.0 )
+        , verbose_( verbose )
+    {}
+
+    General3rdOrderController
+    General3rdOrderController::serializationTestObject()
+    {
+        General3rdOrderController result(1.0, 0.8, true);
+        result.errors_ = {2.0, 3.0};
+
+        return result;
+    }
+
+    double General3rdOrderController::
+    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 )
+        {
+            errors_[i] = errors_[i+1];
+            timeSteps_[i] = timeSteps_[i+1];
+        }
+
+        // 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.)
+        {
+            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_ )
+                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_ )
+                OpmLog::info(fmt::format("Computed step size (pow): {} days", unit::convert::to( newDt, unit::day )));
+            return newDt;
+        }
+    }
+
+    bool General3rdOrderController::operator==(const General3rdOrderController& ctrl) const
+    {
+        return this->tolerance_ == ctrl.tolerance_ &&
+               this->safetyFactor_ == ctrl.safetyFactor_ &&
+               this->errors_ == ctrl.errors_ &&
+               this->timeSteps_ == ctrl.timeSteps_ &&
+               this->verbose_ == ctrl.verbose_;
+    }
+
 } // end namespace Opm

opm/simulators/timestepping/TimeStepControl.hpp

@@ -22,6 +22,7 @@
 #define OPM_TIMESTEPCONTROL_HEADER_INCLUDED
 
 #include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
+#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
 
 #include <string>
 #include <vector>
@@ -32,7 +33,8 @@ namespace Opm
       SimpleIterationCount,
       PID,
       PIDAndIterationCount,
-      HardCodedTimeStep
+      HardCodedTimeStep,
+      General3rdOrder,
     };
 
     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -62,7 +64,7 @@ namespace Opm
         double computeTimeStepSize(const double dt,
                                    const int iterations,
                                    const RelativeChangeInterface& /* relativeChange */,
-                                   const double /*simulationTimeElapsed */ ) const override;
+                                   const AdaptiveSimulatorTimer& /* substepTimer */ ) const override;
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -113,7 +115,7 @@ namespace Opm
         double computeTimeStepSize(const double dt,
                                    const int /* iterations */,
                                    const RelativeChangeInterface& relativeChange,
-                                   const double /*simulationTimeElapsed */ ) const override;
+                                   const AdaptiveSimulatorTimer& /* substepTimer */ ) const override;
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -162,7 +164,7 @@ namespace Opm
         double computeTimeStepSize(const double dt,
                                    const int iterations,
                                    const RelativeChangeInterface& relativeChange,
-                                   const double /*simulationTimeElapsed */ ) const override;
+                                   const AdaptiveSimulatorTimer& /* substepTimer */ ) const override;
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)
@@ -183,6 +185,49 @@ namespace Opm
         const double  minTimeStepBasedOnIterations_;
     };
 
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    ///
+    ///  General 3rd order controller
+    ///
+    ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    class General3rdOrderController : public TimeStepControlInterface
+    {
+    public:
+        static constexpr TimeStepControlType Type = TimeStepControlType::General3rdOrder;
+
+        General3rdOrderController( const double tolerance = 1e-3,
+                                   const double safetyFactor = 0.8,
+                                   const bool verbose = false );
+
+        static General3rdOrderController serializationTestObject();
+
+        double computeTimeStepSize(const double dt,
+                                   const int /*iterations */,
+                                   const RelativeChangeInterface& relativeChange,
+                                   const AdaptiveSimulatorTimer& substepTimer) const override;
+
+        template<class Serializer>
+        void serializeOp(Serializer& serializer)
+        {
+            serializer(tolerance_);
+            serializer(safetyFactor_);
+            serializer(errors_);
+            serializer(timeSteps_);
+            serializer(verbose_);
+        }
+
+        bool operator==(const General3rdOrderController&) const;
+
+
+    protected:
+        const double tolerance_ = 1e-3;
+        const double safetyFactor_ = 0.8;
+        mutable std::vector<double> errors_{};
+        mutable std::vector<double> timeSteps_{};
+        mutable int counterSinceFailure_ = 0;
+        const bool verbose_ = false;
+    };
+
     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ///
     ///  HardcodedTimeStepControl
@@ -208,7 +253,7 @@ namespace Opm
         double computeTimeStepSize(const double dt,
                                    const int /* iterations */,
                                    const RelativeChangeInterface& /*relativeChange */,
-                                   const double simulationTimeElapsed) const override;
+                                   const AdaptiveSimulatorTimer& substepTimer) const override;
 
         template<class Serializer>
         void serializeOp(Serializer& serializer)

opm/simulators/timestepping/TimeStepControlInterface.hpp

@@ -19,6 +19,8 @@
 #ifndef OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
 #define OPM_TIMESTEPCONTROLINTERFACE_HEADER_INCLUDED
 
+#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
+
 
 namespace Opm
 {
@@ -56,7 +58,7 @@ namespace Opm
         /// \param timeError   object to compute || u^n+1 - u^n || / || u^n+1 ||
         ///
         /// \return suggested time step size for the next step
-        virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange , const double simulationTimeElapsed) const = 0;
+        virtual double computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relativeChange , const AdaptiveSimulatorTimer& substepTimer) const = 0;
 
         /// virtual destructor (empty)
         virtual ~TimeStepControlInterface () {}

tests/test_RestartSerialization.cpp

@@ -123,6 +123,7 @@ TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingHardcoded, serializationTestObj
 TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingPID, serializationTestObjectPID)
 TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingPIDIt, serializationTestObjectPIDIt)
 TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeSteppingSimple, serializationTestObjectSimple)
+TEST_FOR_TYPE_NAMED_OBJ(ATS, AdaptiveTimeStepping3rdOrder, serializationTestObject3rdOrder)
 
 namespace Opm { using BPV = BlackOilPrimaryVariables<Properties::TTag::TestTypeTag>; }
 TEST_FOR_TYPE_NAMED(BPV, BlackoilPrimaryVariables)