Merge pull request #2827 from totto82/timestepping

add options for timestepping

opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp

@@ -93,6 +93,14 @@ struct TimeStepControlGrowthRate {
     using type = UndefinedProperty;
 };
 template<class TypeTag, class MyTypeTag>
+struct TimeStepControlDecayDampingFactor {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct TimeStepControlGrowthDampingFactor {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
 struct TimeStepControlFileName {
     using type = UndefinedProperty;
 };
@@ -180,6 +188,16 @@ struct TimeStepControlGrowthRate<TypeTag, TTag::FlowTimeSteppingParameters> {
     static constexpr type value = 1.25;
 };
 template<class TypeTag>
+struct TimeStepControlDecayDampingFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1.0;
+};
+template<class TypeTag>
+struct TimeStepControlGrowthDampingFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1.0/1.2;
+};
+template<class TypeTag>
 struct TimeStepControlFileName<TypeTag, TTag::FlowTimeSteppingParameters> {
     static constexpr auto value = "timesteps";
 };
@@ -296,7 +314,7 @@ namespace Opm {
             EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepAfterEventInDays,
                                  "Time step size of the first time step after an event occurs during the simulation in days");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControl,
-                                 "The algorithm used to determine time-step sizes. valid options are: 'pid' (default), 'pid+iteration', 'pid+newtoniteration', 'iterationcount' and 'hardcoded'");
+                                 "The algorithm used to determine time-step sizes. valid options are: 'pid' (default), 'pid+iteration', 'pid+newtoniteration', 'iterationcount', 'newtoniterationcount' and 'hardcoded'");
             EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlTolerance,
                                  "The tolerance used by the time step size control algorithm");
             EWOMS_REGISTER_PARAM(TypeTag, int, TimeStepControlTargetIterations,
@@ -307,6 +325,10 @@ namespace Opm {
                                  "The decay rate of the time step size of the number of target iterations is exceeded");
             EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthRate,
                                  "The growth rate of the time step size of the number of target iterations is undercut");
+            EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor,
+                                 "The decay rate of the time step decrease when the target iterations is exceeded");
+            EWOMS_REGISTER_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor,
+                                 "The growth rate of the time step increase when the target iterations is undercut");
             EWOMS_REGISTER_PARAM(TypeTag, std::string, TimeStepControlFileName,
                                  "The name of the file which contains the hardcoded time steps sizes");
             EWOMS_REGISTER_PARAM(TypeTag, double, MinTimeStepBeforeShuttingProblematicWellsInDays,
@@ -606,11 +628,15 @@ namespace Opm {
             }
             else if (control == "pid+iteration") {
                 const int iterations =  EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetIterations); // 30
-                timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
+                const double decayDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor); // 1.0
+                const double growthDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor); // 1.0/1.2
+                timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, decayDampingFactor, growthDampingFactor, tol));
             }
             else if (control == "pid+newtoniteration") {
                 const int iterations =  EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
-                timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, tol));
+                const double decayDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayDampingFactor); // 1.0
+                const double growthDampingFactor = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthDampingFactor); // 1.0/1.2
+                timeStepControl_ = TimeStepControlType(new PIDAndIterationCountTimeStepControl(iterations, decayDampingFactor, growthDampingFactor, tol));
                 useNewtonIteration_ = true;
             }
             else if (control == "iterationcount") {
@@ -619,6 +645,13 @@ namespace Opm {
                 const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
                 timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
             }
+            else if (control == "newtoniterationcount") {
+                const int iterations =  EWOMS_GET_PARAM(TypeTag, int, TimeStepControlTargetNewtonIterations); // 8
+                const double decayrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlDecayRate); // 0.75
+                const double growthrate = EWOMS_GET_PARAM(TypeTag, double, TimeStepControlGrowthRate); // 1.25
+                timeStepControl_ = TimeStepControlType(new SimpleIterationCountTimeStepControl(iterations, decayrate, growthrate));
+                useNewtonIteration_ = true;
+            }
             else if (control == "hardcoded") {
                 const std::string filename = EWOMS_GET_PARAM(TypeTag, std::string, TimeStepControlFileName); // "timesteps"
                 timeStepControl_ = TimeStepControlType(new HardcodedTimeStepControl(filename));

opm/simulators/timestepping/TimeStepControl.cpp

@@ -172,26 +172,30 @@ namespace Opm
 
     PIDAndIterationCountTimeStepControl::
     PIDAndIterationCountTimeStepControl( const int target_iterations,
+                                         const double decayDampingFactor,
+                                         const double growthDampingFactor,
                                          const double tol,
                                          const bool verbose)
-        : BaseType( tol, verbose )
+        : PIDTimeStepControl( tol, verbose )
         , target_iterations_( target_iterations )
+        , decayDampingFactor_( decayDampingFactor )
+        , growthDampingFactor_( growthDampingFactor )
     {}
 
     double PIDAndIterationCountTimeStepControl::
     computeTimeStepSize( const double dt, const int iterations, const RelativeChangeInterface& relChange,  const double simulationTimeElapsed ) const
     {
-        double dtEstimatePID = BaseType :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed);
+        double dtEstimatePID = PIDTimeStepControl :: computeTimeStepSize( dt, iterations, relChange, simulationTimeElapsed);
 
         // adjust timesteps based on target iteration
         double dtEstimateIter;
         if (iterations > target_iterations_) {
             double off_target_fraction = double(iterations - target_iterations_) / target_iterations_;
-            dtEstimateIter = dt / (1.0 + off_target_fraction);
+            dtEstimateIter = dt / (1.0 + off_target_fraction * decayDampingFactor_);
         } else {
             double off_target_fraction = double(target_iterations_ - iterations) / target_iterations_;
-            // Be a bit more careful when increasing. The 1.2 factor is from ebos.
+            // Be a bit more careful when increasing.
-            dtEstimateIter = dt * (1.0 + off_target_fraction / 1.2);
+            dtEstimateIter = dt * (1.0 + off_target_fraction * growthDampingFactor_);
         }
 
         return std::min(dtEstimatePID, dtEstimateIter);

opm/simulators/timestepping/TimeStepControl.hpp

@@ -106,6 +106,8 @@ namespace Opm
         ///                   in one time step (default is 1e-3)
         /// \param verbose    if true get some output (default = false)
         PIDAndIterationCountTimeStepControl( const int target_iterations = 20,
+                                             const double decayDampingFactor = 1.0,
+                                             const double growthDampingFactor = 1.0/1.2,
                                              const double tol = 1e-3,
                                              const bool verbose = false);
 
@@ -114,6 +116,8 @@ namespace Opm
 
     protected:
         const int     target_iterations_;
+        const double  decayDampingFactor_;
+        const double  growthDampingFactor_;
     };
 
     ///////////////////////////////////////////////////////////////////////////////////////////////////////////////