diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index aefd590d2526..f410a810e075 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -19,21 +19,25 @@ namespace Kratos { -AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime, - double EndTime, - double StartIncrement, - std::size_t MaxNumOfCycles, - double ReductionFactor, - double IncreaseFactor, - double MaxTimeStepFactor, - std::size_t MinNumOfIterations, - std::size_t MaxNumOfIterations) +AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime, + double EndTime, + double StartIncrement, + std::size_t MaxNumOfCycles, + double ReductionFactor, + double IncreaseFactor, + std::optional mUserMinDeltaTime, + double MaxTimeStepFactor, + std::size_t MinNumOfIterations, + std::size_t MaxNumOfIterations) : TimeIncrementor(), mEndTime(EndTime), - mDeltaTime(std::min(StartIncrement, EndTime - StartTime)), // avoid exceeding the end time + mTimeSpan(EndTime - StartTime), + mDeltaTime(std::min(StartIncrement, mTimeSpan)), // avoid exceeding the end time + mInitialDeltaTime(mDeltaTime), mMaxNumOfCycles(MaxNumOfCycles), mReductionFactor(ReductionFactor), mIncreaseFactor(IncreaseFactor), + mUserMinDeltaTime(std::move(mUserMinDeltaTime)), mMaxDeltaTime(MaxTimeStepFactor * mDeltaTime), mMinNumOfIterations(MinNumOfIterations), mMaxNumOfIterations(MaxNumOfIterations) @@ -73,21 +77,32 @@ 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); - } + constexpr auto delta_time_as_fraction_of_time_span = 0.001; + auto default_min_delta_time = std::min(mInitialDeltaTime, delta_time_as_fraction_of_time_span * mTimeSpan); + + if (rResultantState.Converged()) { + if (rResultantState.num_of_iterations < mMinNumOfIterations) { + mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); + } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) { + mDeltaTime *= mReductionFactor; + } - // Avoid incrementing the time beyond the end time - mDeltaTime = std::min(mDeltaTime, mEndTime - rResultantState.time); + if (rResultantState.time < mEndTime && mEndTime - (rResultantState.time + mDeltaTime) < + mUserMinDeltaTime.value_or(default_min_delta_time)) { + // Up-scaling to reach end_time without small increment + 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 < mUserMinDeltaTime.value_or(default_min_delta_time)) + << "Delta time (" << mDeltaTime << ") is smaller than " + << (mUserMinDeltaTime ? "given" : "default") << " minimum allowable value " + << mUserMinDeltaTime.value_or(default_min_delta_time) << std::endl; } } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h index f99177368849..b1ba9acd8390 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h @@ -15,6 +15,7 @@ #include "includes/kratos_export_api.h" #include "time_incrementor.h" +#include namespace Kratos { @@ -22,15 +23,16 @@ 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::size_t MaxNumOfCycles = 10, + double ReductionFactor = 0.5, + double IncreaseFactor = 2.0, + std::optional mUserMinDeltaTime = std::nullopt, + 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; @@ -38,14 +40,17 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public Tim 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 mTimeSpan; + double mDeltaTime; + double mInitialDeltaTime; + std::size_t mMaxNumOfCycles; + double mReductionFactor; + double mIncreaseFactor; + std::optional mUserMinDeltaTime; + double mMaxDeltaTime; + std::size_t mMinNumOfIterations; + std::size_t mMaxNumOfIterations; }; } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp index 4447f27177fc..1b2d00949615 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp @@ -75,6 +75,15 @@ double GetMaxDeltaTimeFactorFrom(const Parameters& rProjectParameters) return rProjectParameters["solver_settings"]["time_stepping"]["max_delta_time_factor"].GetDouble(); } +std::optional GetMaybeMinDeltaTimeFrom(const Parameters& rProjectParameters) +{ + return rProjectParameters["solver_settings"]["time_stepping"].Has("minimum_allowable_value") + ? std::make_optional(rProjectParameters["solver_settings"]["time_stepping"] + ["minimum_allowable_value"] + .GetDouble()) + : std::nullopt; +} + std::size_t GetMinNumberOfIterationsFrom(const Parameters& rProjectParameters) { return static_cast(rProjectParameters["solver_settings"]["min_iterations"].GetInt()); @@ -325,8 +334,8 @@ std::unique_ptr KratosGeoSettlement::MakeTimeIncrementor(const GetStartTimeFrom(rProjectParameters), GetEndTimeFrom(rProjectParameters), GetTimeIncrementFrom(rProjectParameters), GetMaxNumberOfCyclesFrom(rProjectParameters), GetReductionFactorFrom(rProjectParameters), GetIncreaseFactorFrom(rProjectParameters), - GetMaxDeltaTimeFactorFrom(rProjectParameters), GetMinNumberOfIterationsFrom(rProjectParameters), - GetMaxNumberOfIterationsFrom(rProjectParameters)); + GetMaybeMinDeltaTimeFrom(rProjectParameters), GetMaxDeltaTimeFactorFrom(rProjectParameters), + GetMinNumberOfIterationsFrom(rProjectParameters), GetMaxNumberOfIterationsFrom(rProjectParameters)); } std::shared_ptr KratosGeoSettlement::MakeStrategyWrapper(const Parameters& rProjectParameters, diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index c1485271d23d..fa7126c861c5 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -34,11 +34,13 @@ def __init__(self, model, project_parameters): solver_settings = project_parameters["solver_settings"] self.delta_time = min(solver_settings["time_stepping"]["time_step"].GetDouble(), self.end_time - self.start_time) + self.initial_delta_time = self.delta_time self.reduction_factor = solver_settings["reduction_factor"].GetDouble() self.increase_factor = solver_settings["increase_factor"].GetDouble() 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 = solver_settings["time_stepping"]["minimum_allowable_value"].GetDouble() if solver_settings["time_stepping"].Has("minimum_allowable_value") else None 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,15 @@ 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 _CheckDeltaTimeSize(self): + min_delta_time = self._GetMinDeltaTimeValueOrDefault() + if self.delta_time < min_delta_time: + origin_of_value = "given" if self.min_delta_time is not None else "default" + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), f"The time step {self.delta_time} is smaller than a {origin_of_value} minimum value of {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 @@ -96,12 +106,14 @@ def RunSolutionLoop(self): new_time = t + self.delta_time # avoid very small remaining time steps - small_time_increment = 1.E-3 * self.delta_time - if self.end_time - new_time < small_time_increment: + if self.end_time - new_time < self._GetMinDeltaTimeValueOrDefault(): new_time = self.end_time self.delta_time = new_time - t KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "Up-scaling to reach end_time without small increments: ", self.delta_time) + print("1: ", self.delta_time, t, new_time) + self._CheckDeltaTimeSize() + # start the new step self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP] += 1 self._GetSolver().main_model_part.CloneTimeStep(new_time) @@ -146,6 +158,8 @@ 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 + print("2: ", self.delta_time, t, new_time) + self._CheckDeltaTimeSize() # Reset displacements to the initial KratosMultiphysics.VariableUtils().UpdateCurrentPosition(self._GetSolver().GetComputingModelPart().Nodes, KratosMultiphysics.DISPLACEMENT,1) for node in self._GetSolver().GetComputingModelPart().Nodes: @@ -153,7 +167,7 @@ def RunSolutionLoop(self): 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): @@ -200,6 +214,10 @@ def _GetOrderOfProcessesInitialization(self): def _GetSimulationName(self): return "GeoMechanics Analysis" + def _GetMinDeltaTimeValueOrDefault(self): + delta_time_as_fraction_of_time_span = 0.001 + return self.min_delta_time if self.min_delta_time is not None else min(self.initial_delta_time,delta_time_as_fraction_of_time_span*(self.end_time - self.start_time)) + if __name__ == '__main__': from sys import argv @@ -210,7 +228,7 @@ def _GetSimulationName(self): err_msg += ' "python geomechanics_analysis.py"\n' err_msg += '- With custom parameter file:\n' err_msg += ' "python geomechanics_analysis.py .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] diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_workflows/test_time_incrementor.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_workflows/test_time_incrementor.cpp index 4cae98e80a2d..b9300ec27b2a 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_workflows/test_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_workflows/test_time_incrementor.cpp @@ -9,6 +9,8 @@ // // Main authors: Wijtze Pieter Kikstra // Anne van de Graaf +// Richard Faasse +// Gennady Markelov // #include "custom_workflows/adaptive_time_incrementor.h" @@ -21,15 +23,16 @@ 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::optional MaybeMinDeltaTime{std::make_optional(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) @@ -37,8 +40,8 @@ AdaptiveTimeIncrementor MakeAdaptiveTimeIncrementor(const AdaptiveTimeIncremento return AdaptiveTimeIncrementor{rSettings.StartTime, rSettings.EndTime, rSettings.StartIncrement, rSettings.MaxNumOfCycles, rSettings.ReductionFactor, rSettings.IncreaseFactor, - rSettings.MaxDeltaTimeFactor, rSettings.MinNumOfIterations, - rSettings.MaxNumOfIterations}; + rSettings.MaybeMinDeltaTime, rSettings.MaxDeltaTimeFactor, + rSettings.MinNumOfIterations, rSettings.MaxNumOfIterations}; } } // namespace @@ -427,4 +430,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