[GeoMechanicsApplication] DSettlement non-convergence exit takes a lot longer than before

applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp

@@ -19,9 +19,10 @@
 namespace Kratos
 {
 
-AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double      StartTime,
-                                                 double      EndTime,
-                                                 double      StartIncrement,
+AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime,
+                                                 double EndTime,
+                                                 double StartIncrement,
+                                                 std::pair<std::string, double> MinAllowableDeltaTime,
                                                  std::size_t MaxNumOfCycles,
                                                  double      ReductionFactor,
                                                  double      IncreaseFactor,
@@ -31,6 +32,7 @@ AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double      StartTime,
     : TimeIncrementor(),
       mEndTime(EndTime),
       mDeltaTime(std::min(StartIncrement, EndTime - StartTime)), // avoid exceeding the end time
+      mMinAllowableDeltaTime(std::move(MinAllowableDeltaTime)),
       mMaxNumOfCycles(MaxNumOfCycles),
       mReductionFactor(ReductionFactor),
       mIncreaseFactor(IncreaseFactor),
@@ -73,21 +75,31 @@ double AdaptiveTimeIncrementor::GetIncrement() const { return mDeltaTime; }
 
 void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rResultantState)
 {
-    if (rResultantState.NonConverged() ||
-        (rResultantState.Converged() && (rResultantState.num_of_iterations == mMaxNumOfIterations))) {
-        mDeltaTime *= mReductionFactor;
-    } else if (rResultantState.Converged() && (rResultantState.num_of_iterations < mMinNumOfIterations)) {
-        mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime);
+    if (rResultantState.Converged()) // it is converged also at the beginning of the cycles
+    {
+        if (rResultantState.num_of_iterations < mMinNumOfIterations) {
+            mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime);
+        } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) {
+            mDeltaTime *= mReductionFactor;
+        }
+        if (rResultantState.time < mEndTime) {
+            // Up-scaling to reach end_time without small increments
+            constexpr auto small_increment_factor = 1.0e-3;
+            if (auto small_time_increment = small_increment_factor * mDeltaTime;
+                rResultantState.time + mDeltaTime > mEndTime - small_time_increment) {
+                mDeltaTime = mEndTime - rResultantState.time;
+            }
+        }
     }
 
-    // Avoid incrementing the time beyond the end time
-    mDeltaTime = std::min(mDeltaTime, mEndTime - rResultantState.time);
-
-    // Avoid very small remaining time steps
-    const auto small_time_increment = 1.E-3 * mDeltaTime;
-    if ((mEndTime - (rResultantState.time + mDeltaTime)) < small_time_increment) {
-        mDeltaTime = mEndTime - rResultantState.time;
+    else {
+        // non converged, scale down step and restart
+        mDeltaTime *= mReductionFactor;
     }
+
+    KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second)
+        << "Delta time (" << mDeltaTime << ") is smaller than " << mMinAllowableDeltaTime.first << " minimum allowable value "
+        << mMinAllowableDeltaTime.second << std::endl;
 }
 
 } // namespace Kratos

applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h

@@ -15,37 +15,40 @@
 
 #include "includes/kratos_export_api.h"
 #include "time_incrementor.h"
+#include <string>
 
 namespace Kratos
 {
 
 class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public TimeIncrementor
 {
 public:
-    AdaptiveTimeIncrementor(double      StartTime,
-                            double      EndTime,
-                            double      StartIncrement,
-                            std::size_t MaxNumOfCycles     = 10,
-                            double      ReductionFactor    = 0.5,
-                            double      IncreaseFactor     = 2.0,
-                            double      MaxTimeStepFactor  = 1000.0,
-                            std::size_t MinNumOfIterations = 3,
-                            std::size_t MaxNumOfIterations = 15);
+    AdaptiveTimeIncrementor(double                         StartTime,
+                            double                         EndTime,
+                            double                         StartIncrement,
+                            std::pair<std::string, double> MinAllowableDeltaTime,
+                            std::size_t                    MaxNumOfCycles     = 10,
+                            double                         ReductionFactor    = 0.5,
+                            double                         IncreaseFactor     = 2.0,
+                            double                         MaxTimeStepFactor  = 1000.0,
+                            std::size_t                    MinNumOfIterations = 3,
+                            std::size_t                    MaxNumOfIterations = 15);
 
     [[nodiscard]] bool WantNextStep(const TimeStepEndState& rPreviousState) const override;
     [[nodiscard]] bool WantRetryStep(std::size_t CycleNumber, const TimeStepEndState& rPreviousState) const override;
     [[nodiscard]] double GetIncrement() const override;
     void                 PostTimeStepExecution(const TimeStepEndState& rResultantState) override;
 
 private:
-    double      mEndTime;
-    double      mDeltaTime;
-    std::size_t mMaxNumOfCycles;
-    double      mReductionFactor;
-    double      mIncreaseFactor;
-    double      mMaxDeltaTime;
-    std::size_t mMinNumOfIterations;
-    std::size_t mMaxNumOfIterations;
+    double                         mEndTime;
+    double                         mDeltaTime;
+    std::pair<std::string, double> mMinAllowableDeltaTime;
+    std::size_t                    mMaxNumOfCycles;
+    double                         mReductionFactor;
+    double                         mIncreaseFactor;
+    double                         mMaxDeltaTime;
+    std::size_t                    mMinNumOfIterations;
+    std::size_t                    mMaxNumOfIterations;
 };
 
 } // namespace Kratos

applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp

@@ -75,6 +75,17 @@ double GetMaxDeltaTimeFactorFrom(const Parameters& rProjectParameters)
     return rProjectParameters["solver_settings"]["time_stepping"]["max_delta_time_factor"].GetDouble();
 }
 
+std::pair<std::string, double> GetMinAllowableDeltaTimeFrom(const Parameters& rProjectParameters)
+{
+    if (rProjectParameters["solver_settings"]["time_stepping"].Has("minimum_allowable_value")) {
+        return {"given",
+                rProjectParameters["solver_settings"]["time_stepping"]["minimum_allowable_value"].GetDouble()};
+    } else {
+        constexpr auto minimum_allowable_value = 1e-10;
+        return {"default", minimum_allowable_value};
+    }
+}
+
 std::size_t GetMinNumberOfIterationsFrom(const Parameters& rProjectParameters)
 {
     return static_cast<std::size_t>(rProjectParameters["solver_settings"]["min_iterations"].GetInt());
@@ -323,10 +334,10 @@ std::unique_ptr<TimeIncrementor> KratosGeoSettlement::MakeTimeIncrementor(const
     // For now, we can create adaptive time incrementors only
     return std::make_unique<AdaptiveTimeIncrementor>(
         GetStartTimeFrom(rProjectParameters), GetEndTimeFrom(rProjectParameters),
-        GetTimeIncrementFrom(rProjectParameters), GetMaxNumberOfCyclesFrom(rProjectParameters),
-        GetReductionFactorFrom(rProjectParameters), GetIncreaseFactorFrom(rProjectParameters),
-        GetMaxDeltaTimeFactorFrom(rProjectParameters), GetMinNumberOfIterationsFrom(rProjectParameters),
-        GetMaxNumberOfIterationsFrom(rProjectParameters));
+        GetTimeIncrementFrom(rProjectParameters), GetMinAllowableDeltaTimeFrom(rProjectParameters),
+        GetMaxNumberOfCyclesFrom(rProjectParameters), GetReductionFactorFrom(rProjectParameters),
+        GetIncreaseFactorFrom(rProjectParameters), GetMaxDeltaTimeFactorFrom(rProjectParameters),
+        GetMinNumberOfIterationsFrom(rProjectParameters), GetMaxNumberOfIterationsFrom(rProjectParameters));
 }
 
 std::shared_ptr<StrategyWrapper> KratosGeoSettlement::MakeStrategyWrapper(const Parameters& rProjectParameters,

applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py

@@ -39,6 +39,8 @@ def __init__(self, model, project_parameters):
         self.min_iterations      = solver_settings["min_iterations"].GetInt()
         self.max_delta_time_factor = solver_settings["time_stepping"]["max_delta_time_factor"].GetDouble() if solver_settings["time_stepping"].Has("max_delta_time_factor") else 1000.0
         self.max_delta_time      = self.delta_time * self.max_delta_time_factor
+        self.min_delta_time_set  = solver_settings["time_stepping"].Has("minimum_allowable_value")
+        self.min_delta_time      = solver_settings["time_stepping"]["minimum_allowable_value"].GetDouble() if self.min_delta_time_set else 1e-10
         self.number_cycles       = solver_settings["number_cycles"].GetInt()
         self.max_iterations      = solver_settings["max_iterations"].GetInt()
         self.solution_type       = solver_settings["solution_type"].GetString()
@@ -70,7 +72,16 @@ def FinalizeSolutionStep(self):
         super().FinalizeSolutionStep()
 
         if self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.NL_ITERATION_NUMBER] > self.max_iterations:
-            raise Exception("max_number_of_iterations_exceeded")
+            raise RuntimeError("max_number_of_iterations_exceeded")
+
+    def _check_delta_time_size(self):
+        if self.delta_time < self.min_delta_time:
+            if self.min_delta_time_set:
+                KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "The time step ", self.delta_time, " is smaller than a given minimum value of ", self.min_delta_time)
+            else:
+                KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "The time step ", self.delta_time, " is smaller than a default minimum value of ", self.min_delta_time)
+            KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "Please check settings in Project Parameters and Materials Files.")
+            raise RuntimeError('The time step is too small!')
 
     def RunSolutionLoop(self):
         """This function executes the solution loop of the AnalysisStage
@@ -102,6 +113,8 @@ def RunSolutionLoop(self):
                 self.delta_time = new_time - t
                 KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "Up-scaling to reach end_time without small increments: ", self.delta_time)
 
+            self._check_delta_time_size()
+
             # start the new step
             self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP] += 1
             self._GetSolver().main_model_part.CloneTimeStep(new_time)
@@ -146,14 +159,15 @@ def RunSolutionLoop(self):
                     # scale down step and restart
                     KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "Down-scaling with factor: ", self.reduction_factor)
                     self.delta_time *= self.reduction_factor
+                    self._check_delta_time_size()
                     # Reset displacements to the initial
                     KratosMultiphysics.VariableUtils().UpdateCurrentPosition(self._GetSolver().GetComputingModelPart().Nodes, KratosMultiphysics.DISPLACEMENT,1)
                     for node in self._GetSolver().GetComputingModelPart().Nodes:
                         dold = node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT,1)
                         node.SetSolutionStepValue(KratosMultiphysics.DISPLACEMENT, 0, dold)
 
             if not converged:
-                raise Exception('The maximum number of cycles is reached without convergence!')
+                raise RuntimeError('The maximum number of cycles is reached without convergence!')
 
             if self._GetSolver().settings["reset_displacements"].GetBool():
                 for idx, node in enumerate(self._GetSolver().GetComputingModelPart().Nodes):
@@ -210,7 +224,7 @@ def _GetSimulationName(self):
         err_msg += '    "python geomechanics_analysis.py"\n'
         err_msg += '- With custom parameter file:\n'
         err_msg += '    "python geomechanics_analysis.py <my-parameter-file>.json"\n'
-        raise Exception(err_msg)
+        raise TypeError(err_msg)
 
     if len(argv) == 2: # ProjectParameters is being passed from outside
         parameter_file_name = argv[1]

applications/GeoMechanicsApplication/tests/cpp_tests/custom_workflows/test_time_incrementor.cpp

@@ -9,6 +9,7 @@
 //
 //  Main authors:    Wijtze Pieter Kikstra
 //                   Anne van de Graaf
+//                   Richard Faasse
 //
 
 #include "custom_workflows/adaptive_time_incrementor.h"
@@ -21,24 +22,25 @@ namespace
 {
 
 struct AdaptiveTimeIncrementorSettings {
-    double      StartTime{0.0};
-    double      EndTime{8.0};
-    double      StartIncrement{0.5};
-    std::size_t MaxNumOfCycles{8};
-    double      ReductionFactor{0.5};
-    double      IncreaseFactor{2.0};
-    double      MaxDeltaTimeFactor{1000.0};
-    std::size_t MinNumOfIterations{3};
-    std::size_t MaxNumOfIterations{15};
+    double                         StartTime{0.0};
+    double                         EndTime{8.0};
+    double                         StartIncrement{0.5};
+    std::pair<std::string, double> MinAllowableDeltaTime{"given", 1e-06};
+    std::size_t                    MaxNumOfCycles{8};
+    double                         ReductionFactor{0.5};
+    double                         IncreaseFactor{2.0};
+    double                         MaxDeltaTimeFactor{1000.0};
+    std::size_t                    MinNumOfIterations{3};
+    std::size_t                    MaxNumOfIterations{15};
 };
 
 AdaptiveTimeIncrementor MakeAdaptiveTimeIncrementor(const AdaptiveTimeIncrementorSettings& rSettings)
 {
     return AdaptiveTimeIncrementor{rSettings.StartTime,          rSettings.EndTime,
-                                   rSettings.StartIncrement,     rSettings.MaxNumOfCycles,
-                                   rSettings.ReductionFactor,    rSettings.IncreaseFactor,
-                                   rSettings.MaxDeltaTimeFactor, rSettings.MinNumOfIterations,
-                                   rSettings.MaxNumOfIterations};
+                                   rSettings.StartIncrement,     rSettings.MinAllowableDeltaTime,
+                                   rSettings.MaxNumOfCycles,     rSettings.ReductionFactor,
+                                   rSettings.IncreaseFactor,     rSettings.MaxDeltaTimeFactor,
+                                   rSettings.MinNumOfIterations, rSettings.MaxNumOfIterations};
 }
 
 } // namespace
@@ -427,4 +429,38 @@ KRATOS_TEST_CASE_IN_SUITE(ScaleIncrementToAvoidExtraSmallTimeStep, KratosGeoMech
     KRATOS_EXPECT_DOUBLE_EQ(8.0, time_incrementor.GetIncrement());
 }
 
+KRATOS_TEST_CASE_IN_SUITE(ThrowExceptionWhenDeltaTimeSmallerThanTheLimit, KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    AdaptiveTimeIncrementorSettings settings; // with EndTime = 8.0
+    settings.StartIncrement          = 8.0;   // We jump to the end time right away
+    auto time_incrementor            = MakeAdaptiveTimeIncrementor(settings);
+    auto previous_state              = TimeStepEndState{};
+    previous_state.convergence_state = TimeStepEndState::ConvergenceState::converged;
+    previous_state.time              = 7.9999999; // to have a zero time step
+
+    KRATOS_EXPECT_EXCEPTION_IS_THROWN(
+        time_incrementor.PostTimeStepExecution(previous_state),
+        "Delta time (1e-07) is smaller than given minimum allowable value 1e-06");
+
+    previous_state.convergence_state = TimeStepEndState::ConvergenceState::non_converged;
+
+    KRATOS_EXPECT_EXCEPTION_IS_THROWN(
+        time_incrementor.PostTimeStepExecution(previous_state),
+        "Delta time (5e-08) is smaller than given minimum allowable value 1e-06");
+}
+
+KRATOS_TEST_CASE_IN_SUITE(HalfTimeStepAtNonConverged, KratosGeoMechanicsFastSuiteWithoutKernel)
+{
+    AdaptiveTimeIncrementorSettings settings; // with EndTime = 8.0
+    settings.StartIncrement          = 8.0;   // We jump to the end time right away
+    auto time_incrementor            = MakeAdaptiveTimeIncrementor(settings);
+    auto previous_state              = TimeStepEndState{};
+    previous_state.convergence_state = TimeStepEndState::ConvergenceState::non_converged;
+    previous_state.time              = 8.0;
+
+    time_incrementor.PostTimeStepExecution(previous_state);
+    // The increment should be halved, since the step didn't converge
+    KRATOS_EXPECT_DOUBLE_EQ(4.0, time_incrementor.GetIncrement());
+}
+
 } // namespace Kratos::Testing