From 3df6c93bdea738b5171b2905e9460c46ebb6b170 Mon Sep 17 00:00:00 2001 From: markelov Date: Tue, 25 Feb 2025 13:36:16 +0100 Subject: [PATCH 01/16] added unit test --- .../custom_workflows/test_time_incrementor.cpp | 15 +++++++++++++++ 1 file changed, 15 insertions(+) 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..a5bf74088f52 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,7 @@ // // Main authors: Wijtze Pieter Kikstra // Anne van de Graaf +// Richard Faasse // #include "custom_workflows/adaptive_time_incrementor.h" @@ -427,4 +428,18 @@ KRATOS_TEST_CASE_IN_SUITE(ScaleIncrementToAvoidExtraSmallTimeStep, KratosGeoMech KRATOS_EXPECT_DOUBLE_EQ(8.0, time_incrementor.GetIncrement()); } +KRATOS_TEST_CASE_IN_SUITE(IncrementIsCorrectWhenNonConvergedStepIsSameAsEndTime, 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; // The jumped step didn't converge + previous_state.time = 8.0; // The non-converged step is the same as the end time + + 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 From af922ed3b9c3b9754c752be355344840aab7e77d Mon Sep 17 00:00:00 2001 From: markelov Date: Fri, 28 Feb 2025 10:59:57 +0100 Subject: [PATCH 02/16] added minimum allowable value for the time step --- .../adaptive_time_incrementor.cpp | 13 ++++-- .../adaptive_time_incrementor.h | 38 +++++++++------- .../custom_workflows/dgeosettlement.cpp | 21 +++++++-- .../python_scripts/geomechanics_analysis.py | 20 +++++++-- .../test_time_incrementor.cpp | 44 +++++++++---------- 5 files changed, 87 insertions(+), 49 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index aefd590d2526..4e7b361f4677 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/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 MinAllowableDeltaTime, std::size_t MaxNumOfCycles, double ReductionFactor, double IncreaseFactor, @@ -29,8 +30,10 @@ AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime, std::size_t MinNumOfIterations, std::size_t MaxNumOfIterations) : TimeIncrementor(), + mTimeSpan(EndTime - StartTime), mEndTime(EndTime), mDeltaTime(std::min(StartIncrement, EndTime - StartTime)), // avoid exceeding the end time + mMinAllowableDeltaTime(std::move(MinAllowableDeltaTime)), mMaxNumOfCycles(MaxNumOfCycles), mReductionFactor(ReductionFactor), mIncreaseFactor(IncreaseFactor), @@ -88,6 +91,10 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes if ((mEndTime - (rResultantState.time + mDeltaTime)) < small_time_increment) { mDeltaTime = mEndTime - rResultantState.time; } + + KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) + << "Delta time (" << mDeltaTime << ") is smaller than minimum allowable value " + << mMinAllowableDeltaTime.second << 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..e39f9e52a202 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::pair 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; @@ -38,14 +40,16 @@ 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 mTimeSpan; + double mEndTime; + double mDeltaTime; + std::pair mMinAllowableDeltaTime; + std::size_t mMaxNumOfCycles; + double mReductionFactor; + double mIncreaseFactor; + 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..abe2f00da103 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp @@ -75,6 +75,19 @@ double GetMaxDeltaTimeFactorFrom(const Parameters& rProjectParameters) return rProjectParameters["solver_settings"]["time_stepping"]["max_delta_time_factor"].GetDouble(); } +std::pair GetMinAllowableDeltaTimeFrom(const Parameters& rProjectParameters) +{ + if (rProjectParameters["solver_settings"]["time_stepping"].Has("minimum_allowable_value")) { + return {"set", + rProjectParameters["solver_settings"]["time_stepping"]["minimum_allowable_value"].GetDouble()}; + } else { + double minimum_allowable_value = (rProjectParameters["solver_settings"]["end_time"].GetDouble() - + rProjectParameters["solver_settings"]["start_time"].GetDouble()) / + 100000.0; + return {"default", minimum_allowable_value}; + } +} + std::size_t GetMinNumberOfIterationsFrom(const Parameters& rProjectParameters) { return static_cast(rProjectParameters["solver_settings"]["min_iterations"].GetInt()); @@ -323,10 +336,10 @@ std::unique_ptr KratosGeoSettlement::MakeTimeIncrementor(const // For now, we can create adaptive time incrementors only return std::make_unique( 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 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..d3102cd3d351 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/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 (self.end_time - self.start_time)/100000.0 self.number_cycles = solver_settings["number_cycles"].GetInt() self.max_iterations = solver_settings["max_iterations"].GetInt() self.solution_type = solver_settings["solution_type"].GetString() @@ -72,6 +74,15 @@ def FinalizeSolutionStep(self): if self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.NL_ITERATION_NUMBER] > self.max_iterations: raise Exception("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 the case settings.") + raise Exception('The time step is too small!') + def RunSolutionLoop(self): """This function executes the solution loop of the AnalysisStage It can be overridden by derived classes @@ -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,6 +159,7 @@ 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: @@ -180,9 +194,9 @@ def ResetIfHasNodalSolutionStepVariable(self, variable): self._GetSolver().GetComputingModelPart().Nodes, variable, zero_vector, 1) def PrintAnalysisStageProgressInformation(self): - KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "STEP : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP]) - KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "DELTA_TIME: ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.DELTA_TIME]) - KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "TIME : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME]) + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "STEP : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP]) + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "DELTA_TIME : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.DELTA_TIME]) + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "TIME : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME]) def _CreateSolver(self): return geomechanics_solvers_wrapper.CreateSolver(self.model, self.project_parameters) 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 a5bf74088f52..4c49154e241f 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 @@ -22,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 MinAllowableDeltaTime{"set", 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 @@ -428,18 +429,17 @@ KRATOS_TEST_CASE_IN_SUITE(ScaleIncrementToAvoidExtraSmallTimeStep, KratosGeoMech KRATOS_EXPECT_DOUBLE_EQ(8.0, time_incrementor.GetIncrement()); } -KRATOS_TEST_CASE_IN_SUITE(IncrementIsCorrectWhenNonConvergedStepIsSameAsEndTime, KratosGeoMechanicsFastSuiteWithoutKernel) +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::non_converged; // The jumped step didn't converge - previous_state.time = 8.0; // The non-converged step is the same as the end time + settings.StartIncrement = 8.0; // We jump to the end time right away + auto time_incrementor = MakeAdaptiveTimeIncrementor(settings); + auto previous_state = TimeStepEndState{}; + previous_state.time = 8.0; // to have a zero time step - 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()); + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + time_incrementor.PostTimeStepExecution(previous_state), + "Delta time (0) is smaller than minimum allowable value 1e-06"); } } // namespace Kratos::Testing From 800d26023b6d31e8a8281539465b3e9774e77c90 Mon Sep 17 00:00:00 2001 From: markelov Date: Fri, 28 Feb 2025 15:12:13 +0100 Subject: [PATCH 03/16] fixed default minimum allowable time step --- .../custom_workflows/dgeosettlement.cpp | 4 +--- .../python_scripts/geomechanics_analysis.py | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp index abe2f00da103..427821c1b889 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp @@ -81,9 +81,7 @@ std::pair GetMinAllowableDeltaTimeFrom(const Parameters& rP return {"set", rProjectParameters["solver_settings"]["time_stepping"]["minimum_allowable_value"].GetDouble()}; } else { - double minimum_allowable_value = (rProjectParameters["solver_settings"]["end_time"].GetDouble() - - rProjectParameters["solver_settings"]["start_time"].GetDouble()) / - 100000.0; + constexpr auto minimum_allowable_value = 1e-10; return {"default", minimum_allowable_value}; } } diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index d3102cd3d351..03bed2d9034c 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -40,7 +40,7 @@ def __init__(self, model, project_parameters): 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 (self.end_time - self.start_time)/100000.0 + 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() From d3009b185d12be801fbb1325326d974eb04d10af Mon Sep 17 00:00:00 2001 From: markelov Date: Fri, 28 Feb 2025 15:41:35 +0100 Subject: [PATCH 04/16] removed unused mTimeSpan --- .../custom_workflows/adaptive_time_incrementor.cpp | 1 - .../custom_workflows/adaptive_time_incrementor.h | 1 - 2 files changed, 2 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index 4e7b361f4677..ae7914e5ecf0 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -30,7 +30,6 @@ AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime, std::size_t MinNumOfIterations, std::size_t MaxNumOfIterations) : TimeIncrementor(), - mTimeSpan(EndTime - StartTime), mEndTime(EndTime), mDeltaTime(std::min(StartIncrement, EndTime - StartTime)), // avoid exceeding the end time mMinAllowableDeltaTime(std::move(MinAllowableDeltaTime)), diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h index e39f9e52a202..fff2645d57e0 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h @@ -40,7 +40,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public Tim void PostTimeStepExecution(const TimeStepEndState& rResultantState) override; private: - double mTimeSpan; double mEndTime; double mDeltaTime; std::pair mMinAllowableDeltaTime; From 0e6780482dd11822d09a8d555ce7dc31fd4c6174 Mon Sep 17 00:00:00 2001 From: markelov Date: Tue, 4 Mar 2025 13:40:57 +0100 Subject: [PATCH 05/16] changed time step function --- .../adaptive_time_incrementor.cpp | 40 +++++++++++-------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index ae7914e5ecf0..c56206238aa6 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -75,25 +75,33 @@ double AdaptiveTimeIncrementor::GetIncrement() const { return mDeltaTime; } void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rResultantState) { - if (rResultantState.NonConverged() || - (rResultantState.Converged() && (rResultantState.num_of_iterations == mMaxNumOfIterations))) { + // mDeltaTime = std::min(mDeltaTime, mMaxDeltaTime); + auto conv = rResultantState.Converged(); + auto nonconv = rResultantState.NonConverged(); + if (rResultantState.Converged()) // it is converged also at the beginning of the cycles + { + // # scale next step if desired + if (rResultantState.time < mEndTime) { + if (rResultantState.num_of_iterations < mMinNumOfIterations) { + // #scale up next step + mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); + constexpr auto fraction_to_avoid_too_small_step = 1.0e-3; + if (rResultantState.time + mDeltaTime > mEndTime - fraction_to_avoid_too_small_step * mDeltaTime) { + mDeltaTime = mEndTime - rResultantState.time; + } + } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) { + // # converged, but max_iterations reached, scale down next step + mDeltaTime *= mReductionFactor; + } + } + } else { + // # scale down step and restart mDeltaTime *= mReductionFactor; - } else if (rResultantState.Converged() && (rResultantState.num_of_iterations < mMinNumOfIterations)) { - mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); - } - - // 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; + KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) + << "Delta time (" << mDeltaTime << ") is smaller than minimum allowable value " + << mMinAllowableDeltaTime.second << std::endl; } - - KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) - << "Delta time (" << mDeltaTime << ") is smaller than minimum allowable value " - << mMinAllowableDeltaTime.second << std::endl; } } // namespace Kratos From a176a87b7353c71470d404e32345925f0ddf5e80 Mon Sep 17 00:00:00 2001 From: markelov Date: Tue, 4 Mar 2025 13:41:17 +0100 Subject: [PATCH 06/16] fixed code smell --- .../python_scripts/geomechanics_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index 03bed2d9034c..14d6d5b927a5 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -167,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): From 9bd1ee33ebb1284a6135bf34497bdaf181dbecae Mon Sep 17 00:00:00 2001 From: markelov Date: Tue, 4 Mar 2025 20:09:53 +0100 Subject: [PATCH 07/16] removed unused variables --- .../custom_workflows/adaptive_time_incrementor.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index c56206238aa6..dceee89724d1 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -75,9 +75,6 @@ double AdaptiveTimeIncrementor::GetIncrement() const { return mDeltaTime; } void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rResultantState) { - // mDeltaTime = std::min(mDeltaTime, mMaxDeltaTime); - auto conv = rResultantState.Converged(); - auto nonconv = rResultantState.NonConverged(); if (rResultantState.Converged()) // it is converged also at the beginning of the cycles { // # scale next step if desired From 9a936a0dc9a568db6f56c78dbea2fae019e38623 Mon Sep 17 00:00:00 2001 From: markelov Date: Tue, 4 Mar 2025 21:05:49 +0100 Subject: [PATCH 08/16] added unit test --- .../adaptive_time_incrementor.cpp | 24 +++++++++---------- .../test_time_incrementor.cpp | 23 ++++++++++++++---- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index dceee89724d1..6894a7c18bc0 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -77,28 +77,28 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes { if (rResultantState.Converged()) // it is converged also at the beginning of the cycles { - // # scale next step if desired + // scale next step if desired if (rResultantState.time < mEndTime) { if (rResultantState.num_of_iterations < mMinNumOfIterations) { - // #scale up next step + // scale up next step mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); - constexpr auto fraction_to_avoid_too_small_step = 1.0e-3; - if (rResultantState.time + mDeltaTime > mEndTime - fraction_to_avoid_too_small_step * mDeltaTime) { - mDeltaTime = mEndTime - rResultantState.time; - } } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) { - // # converged, but max_iterations reached, scale down next step + // converged, but max_iterations reached, scale down next step mDeltaTime *= mReductionFactor; } } + constexpr auto fraction_to_avoid_too_small_step = 1.0e-3; + if (rResultantState.time + mDeltaTime > mEndTime - fraction_to_avoid_too_small_step * mDeltaTime) { + mDeltaTime = mEndTime - rResultantState.time; + } } else { - // # scale down step and restart + // non converged, scale down step and restart mDeltaTime *= mReductionFactor; - - KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) - << "Delta time (" << mDeltaTime << ") is smaller than minimum allowable value " - << mMinAllowableDeltaTime.second << std::endl; } + + KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) + << "Delta time (" << mDeltaTime << ") is smaller than minimum allowable value " + << mMinAllowableDeltaTime.second << std::endl; } } // namespace Kratos 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 4c49154e241f..3302c3cebe85 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 @@ -432,14 +432,29 @@ KRATOS_TEST_CASE_IN_SUITE(ScaleIncrementToAvoidExtraSmallTimeStep, KratosGeoMech 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.time = 8.0; // to have a zero time step + 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 = 8.0; // to have a zero time step KRATOS_EXPECT_EXCEPTION_IS_THROWN( time_incrementor.PostTimeStepExecution(previous_state), "Delta time (0) is smaller than 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 From 0bca4e1d33fa76913fcd21a7afec931a5b2dfb5a Mon Sep 17 00:00:00 2001 From: markelov Date: Wed, 5 Mar 2025 13:11:39 +0100 Subject: [PATCH 09/16] replaced Exception with other exception types --- .../python_scripts/geomechanics_analysis.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index 14d6d5b927a5..3ed6adc56c74 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -72,7 +72,7 @@ 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: @@ -81,7 +81,7 @@ def _check_delta_time_size(self): 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 the case settings.") - raise Exception('The time step is too small!') + raise RuntimeError('The time step is too small!') def RunSolutionLoop(self): """This function executes the solution loop of the AnalysisStage @@ -224,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 .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] From 32f419898c8bdf1f03d26bdb4081b324f2e013d4 Mon Sep 17 00:00:00 2001 From: markelov Date: Wed, 5 Mar 2025 13:47:40 +0100 Subject: [PATCH 10/16] avoid small increment for non-finished computation --- .../adaptive_time_incrementor.cpp | 23 +++++++++---------- .../test_time_incrementor.cpp | 4 ++-- 2 files changed, 13 insertions(+), 14 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index 6894a7c18bc0..6f3328e95494 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -77,20 +77,19 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes { if (rResultantState.Converged()) // it is converged also at the beginning of the cycles { - // scale next step if desired - if (rResultantState.time < mEndTime) { - if (rResultantState.num_of_iterations < mMinNumOfIterations) { - // scale up next step - mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); - } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) { - // converged, but max_iterations reached, scale down next step - mDeltaTime *= mReductionFactor; - } + if (rResultantState.num_of_iterations < mMinNumOfIterations) { + mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); + } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) { + mDeltaTime *= mReductionFactor; } - constexpr auto fraction_to_avoid_too_small_step = 1.0e-3; - if (rResultantState.time + mDeltaTime > mEndTime - fraction_to_avoid_too_small_step * mDeltaTime) { - mDeltaTime = mEndTime - rResultantState.time; + if (rResultantState.time < mEndTime) { + // Up-scaling to reach end_time without small increments + constexpr auto small_increment = 1.0e-3; + mDeltaTime = (rResultantState.time + mDeltaTime > mEndTime - small_increment * mDeltaTime) + ? mDeltaTime = mEndTime - rResultantState.time + : mDeltaTime; } + } else { // non converged, scale down step and restart mDeltaTime *= mReductionFactor; 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 3302c3cebe85..22e49fdd9b7b 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 @@ -436,11 +436,11 @@ KRATOS_TEST_CASE_IN_SUITE(ThrowExceptionWhenDeltaTimeSmallerThanTheLimit, Kratos auto time_incrementor = MakeAdaptiveTimeIncrementor(settings); auto previous_state = TimeStepEndState{}; previous_state.convergence_state = TimeStepEndState::ConvergenceState::converged; - previous_state.time = 8.0; // to have a zero time step + previous_state.time = 7.9999999; // to have a zero time step KRATOS_EXPECT_EXCEPTION_IS_THROWN( time_incrementor.PostTimeStepExecution(previous_state), - "Delta time (0) is smaller than minimum allowable value 1e-06"); + "Delta time (1e-07) is smaller than minimum allowable value 1e-06"); } KRATOS_TEST_CASE_IN_SUITE(HalfTimeStepAtNonConverged, KratosGeoMechanicsFastSuiteWithoutKernel) From 18e69e1b28acbc8fb08adfdab5b361dc224c2c4b Mon Sep 17 00:00:00 2001 From: markelov Date: Wed, 5 Mar 2025 14:24:12 +0100 Subject: [PATCH 11/16] fixed gcc compilation error --- .../custom_workflows/adaptive_time_incrementor.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index 6f3328e95494..72222063c1e1 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -85,12 +85,13 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes if (rResultantState.time < mEndTime) { // Up-scaling to reach end_time without small increments constexpr auto small_increment = 1.0e-3; - mDeltaTime = (rResultantState.time + mDeltaTime > mEndTime - small_increment * mDeltaTime) - ? mDeltaTime = mEndTime - rResultantState.time - : mDeltaTime; + if (rResultantState.time + mDeltaTime > mEndTime - small_increment * mDeltaTime) { + mDeltaTime = mEndTime - rResultantState.time; + } } + } - } else { + else { // non converged, scale down step and restart mDeltaTime *= mReductionFactor; } From 449ac959f58e48a14936664ffaebce716f7498dd Mon Sep 17 00:00:00 2001 From: markelov Date: Wed, 5 Mar 2025 16:29:55 +0100 Subject: [PATCH 12/16] extended unit test for limit --- .../cpp_tests/custom_workflows/test_time_incrementor.cpp | 6 ++++++ 1 file changed, 6 insertions(+) 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 22e49fdd9b7b..defc8a33e72c 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 @@ -441,6 +441,12 @@ KRATOS_TEST_CASE_IN_SUITE(ThrowExceptionWhenDeltaTimeSmallerThanTheLimit, Kratos KRATOS_EXPECT_EXCEPTION_IS_THROWN( time_incrementor.PostTimeStepExecution(previous_state), "Delta time (1e-07) is smaller than 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 minimum allowable value 1e-06"); } KRATOS_TEST_CASE_IN_SUITE(HalfTimeStepAtNonConverged, KratosGeoMechanicsFastSuiteWithoutKernel) From b6d6e6c9dca4f7cc407d447fa94690a4644a0d4d Mon Sep 17 00:00:00 2001 From: markelov Date: Thu, 6 Mar 2025 12:11:06 +0100 Subject: [PATCH 13/16] response to WPK comments --- .../custom_workflows/adaptive_time_incrementor.cpp | 7 ++++--- .../custom_workflows/dgeosettlement.cpp | 2 +- .../python_scripts/geomechanics_analysis.py | 8 ++++---- .../cpp_tests/custom_workflows/test_time_incrementor.cpp | 8 ++++---- 4 files changed, 13 insertions(+), 12 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index 72222063c1e1..ebcc8fe1f8c6 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -84,8 +84,9 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes } if (rResultantState.time < mEndTime) { // Up-scaling to reach end_time without small increments - constexpr auto small_increment = 1.0e-3; - if (rResultantState.time + mDeltaTime > mEndTime - small_increment * mDeltaTime) { + 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; } } @@ -97,7 +98,7 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes } KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) - << "Delta time (" << mDeltaTime << ") is smaller than minimum allowable value " + << "Delta time (" << mDeltaTime << ") is smaller than " << mMinAllowableDeltaTime.first << " minimum allowable value " << mMinAllowableDeltaTime.second << std::endl; } diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp index 427821c1b889..180f93857393 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp @@ -78,7 +78,7 @@ double GetMaxDeltaTimeFactorFrom(const Parameters& rProjectParameters) std::pair GetMinAllowableDeltaTimeFrom(const Parameters& rProjectParameters) { if (rProjectParameters["solver_settings"]["time_stepping"].Has("minimum_allowable_value")) { - return {"set", + return {"given", rProjectParameters["solver_settings"]["time_stepping"]["minimum_allowable_value"].GetDouble()}; } else { constexpr auto minimum_allowable_value = 1e-10; diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index 3ed6adc56c74..fd8b66b4a9be 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -80,7 +80,7 @@ def _check_delta_time_size(self): 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 the case settings.") + 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): @@ -194,9 +194,9 @@ def ResetIfHasNodalSolutionStepVariable(self, variable): self._GetSolver().GetComputingModelPart().Nodes, variable, zero_vector, 1) def PrintAnalysisStageProgressInformation(self): - KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "STEP : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP]) - KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "DELTA_TIME : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.DELTA_TIME]) - KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "TIME : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME]) + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "STEP : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP]) + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "DELTA_TIME: ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.DELTA_TIME]) + KratosMultiphysics.Logger.PrintInfo(self._GetSimulationName(), "TIME : ", self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.TIME]) def _CreateSolver(self): return geomechanics_solvers_wrapper.CreateSolver(self.model, self.project_parameters) 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 defc8a33e72c..fb526827f90b 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 @@ -25,7 +25,7 @@ struct AdaptiveTimeIncrementorSettings { double StartTime{0.0}; double EndTime{8.0}; double StartIncrement{0.5}; - std::pair MinAllowableDeltaTime{"set", 1e-06}; + std::pair MinAllowableDeltaTime{"given", 1e-06}; std::size_t MaxNumOfCycles{8}; double ReductionFactor{0.5}; double IncreaseFactor{2.0}; @@ -440,13 +440,13 @@ KRATOS_TEST_CASE_IN_SUITE(ThrowExceptionWhenDeltaTimeSmallerThanTheLimit, Kratos KRATOS_EXPECT_EXCEPTION_IS_THROWN( time_incrementor.PostTimeStepExecution(previous_state), - "Delta time (1e-07) is smaller than minimum allowable value 1e-06"); + "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 minimum allowable value 1e-06"); + 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) From 407d2cfc4eeb0c9fa0d04adda39106bb8742c7c5 Mon Sep 17 00:00:00 2001 From: markelov Date: Thu, 6 Mar 2025 23:30:52 +0100 Subject: [PATCH 14/16] response to Anne's comments --- .../adaptive_time_incrementor.cpp | 43 ++++++++++--------- .../adaptive_time_incrementor.h | 38 ++++++++-------- .../custom_workflows/dgeosettlement.cpp | 20 ++++----- .../python_scripts/geomechanics_analysis.py | 26 +++++------ .../test_time_incrementor.cpp | 26 +++++------ 5 files changed, 76 insertions(+), 77 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index ebcc8fe1f8c6..7d1715952888 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -19,23 +19,23 @@ namespace Kratos { -AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime, - double EndTime, - double StartIncrement, - std::pair MinAllowableDeltaTime, - 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 MaybeMinDeltaTime, + double MaxTimeStepFactor, + std::size_t MinNumOfIterations, + std::size_t MaxNumOfIterations) : TimeIncrementor(), mEndTime(EndTime), mDeltaTime(std::min(StartIncrement, EndTime - StartTime)), // avoid exceeding the end time - mMinAllowableDeltaTime(std::move(MinAllowableDeltaTime)), mMaxNumOfCycles(MaxNumOfCycles), mReductionFactor(ReductionFactor), mIncreaseFactor(IncreaseFactor), + mMaybeMinDeltaTime(std::move(MaybeMinDeltaTime)), mMaxDeltaTime(MaxTimeStepFactor * mDeltaTime), mMinNumOfIterations(MinNumOfIterations), mMaxNumOfIterations(MaxNumOfIterations) @@ -75,6 +75,8 @@ double AdaptiveTimeIncrementor::GetIncrement() const { return mDeltaTime; } void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rResultantState) { + constexpr auto default_min_delta_time = 1.0e-10; + if (rResultantState.Converged()) // it is converged also at the beginning of the cycles { if (rResultantState.num_of_iterations < mMinNumOfIterations) { @@ -82,13 +84,11 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes } 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; - } + + if (rResultantState.time < mEndTime && mEndTime - (rResultantState.time + mDeltaTime) < + mMaybeMinDeltaTime.value_or(default_min_delta_time)) { + // Up-scaling to reach end_time without small increment + mDeltaTime = mEndTime - rResultantState.time; } } @@ -97,9 +97,10 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes mDeltaTime *= mReductionFactor; } - KRATOS_ERROR_IF(mDeltaTime < mMinAllowableDeltaTime.second) - << "Delta time (" << mDeltaTime << ") is smaller than " << mMinAllowableDeltaTime.first << " minimum allowable value " - << mMinAllowableDeltaTime.second << std::endl; + KRATOS_ERROR_IF(mDeltaTime < mMaybeMinDeltaTime.value_or(default_min_delta_time)) + << "Delta time (" << mDeltaTime << ") is smaller than " + << (mMaybeMinDeltaTime ? "given" : "default") << " minimum allowable value " + << mMaybeMinDeltaTime.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 fff2645d57e0..911a34ebd991 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h @@ -23,16 +23,16 @@ namespace Kratos class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public TimeIncrementor { public: - AdaptiveTimeIncrementor(double StartTime, - double EndTime, - double StartIncrement, - std::pair 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); + AdaptiveTimeIncrementor(double StartTime, + double EndTime, + double StartIncrement, + std::size_t MaxNumOfCycles = 10, + double ReductionFactor = 0.5, + double IncreaseFactor = 2.0, + std::optional MaybeMinDeltaTime = 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; @@ -40,15 +40,15 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public Tim void PostTimeStepExecution(const TimeStepEndState& rResultantState) override; private: - double mEndTime; - double mDeltaTime; - std::pair mMinAllowableDeltaTime; - std::size_t mMaxNumOfCycles; - double mReductionFactor; - double mIncreaseFactor; - double mMaxDeltaTime; - std::size_t mMinNumOfIterations; - std::size_t mMaxNumOfIterations; + double mEndTime; + double mDeltaTime; + std::size_t mMaxNumOfCycles; + double mReductionFactor; + double mIncreaseFactor; + std::optional mMaybeMinDeltaTime; + 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 180f93857393..1b2d00949615 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp @@ -75,15 +75,13 @@ double GetMaxDeltaTimeFactorFrom(const Parameters& rProjectParameters) return rProjectParameters["solver_settings"]["time_stepping"]["max_delta_time_factor"].GetDouble(); } -std::pair GetMinAllowableDeltaTimeFrom(const Parameters& rProjectParameters) +std::optional GetMaybeMinDeltaTimeFrom(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}; - } + 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) @@ -334,9 +332,9 @@ std::unique_ptr KratosGeoSettlement::MakeTimeIncrementor(const // For now, we can create adaptive time incrementors only return std::make_unique( GetStartTimeFrom(rProjectParameters), GetEndTimeFrom(rProjectParameters), - GetTimeIncrementFrom(rProjectParameters), GetMinAllowableDeltaTimeFrom(rProjectParameters), - GetMaxNumberOfCyclesFrom(rProjectParameters), GetReductionFactorFrom(rProjectParameters), - GetIncreaseFactorFrom(rProjectParameters), GetMaxDeltaTimeFactorFrom(rProjectParameters), + GetTimeIncrementFrom(rProjectParameters), GetMaxNumberOfCyclesFrom(rProjectParameters), + GetReductionFactorFrom(rProjectParameters), GetIncreaseFactorFrom(rProjectParameters), + GetMaybeMinDeltaTimeFrom(rProjectParameters), GetMaxDeltaTimeFactorFrom(rProjectParameters), GetMinNumberOfIterationsFrom(rProjectParameters), GetMaxNumberOfIterationsFrom(rProjectParameters)); } diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index fd8b66b4a9be..2778b822539f 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -39,8 +39,7 @@ 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.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() @@ -74,13 +73,12 @@ def FinalizeSolutionStep(self): if self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.NL_ITERATION_NUMBER] > self.max_iterations: 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.") + 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): @@ -107,13 +105,12 @@ 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) - self._check_delta_time_size() + self._CheckDeltaTimeSize() # start the new step self._GetSolver().GetComputingModelPart().ProcessInfo[KratosMultiphysics.STEP] += 1 @@ -159,7 +156,7 @@ 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() + self._CheckDeltaTimeSize() # Reset displacements to the initial KratosMultiphysics.VariableUtils().UpdateCurrentPosition(self._GetSolver().GetComputingModelPart().Nodes, KratosMultiphysics.DISPLACEMENT,1) for node in self._GetSolver().GetComputingModelPart().Nodes: @@ -214,6 +211,9 @@ def _GetOrderOfProcessesInitialization(self): def _GetSimulationName(self): return "GeoMechanics Analysis" + def _GetMinDeltaTimeValueOrDefault(self): + return self.min_delta_time if self.min_delta_time is not None else 1.0e-10 + if __name__ == '__main__': from sys import argv 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 fb526827f90b..4d5bfe448ec3 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 @@ -22,24 +22,24 @@ namespace { struct AdaptiveTimeIncrementorSettings { - double StartTime{0.0}; - double EndTime{8.0}; - double StartIncrement{0.5}; - std::pair 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}; + 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) { return AdaptiveTimeIncrementor{rSettings.StartTime, rSettings.EndTime, - rSettings.StartIncrement, rSettings.MinAllowableDeltaTime, - rSettings.MaxNumOfCycles, rSettings.ReductionFactor, - rSettings.IncreaseFactor, rSettings.MaxDeltaTimeFactor, + rSettings.StartIncrement, rSettings.MaxNumOfCycles, + rSettings.ReductionFactor, rSettings.IncreaseFactor, + rSettings.MaybeMinDeltaTime, rSettings.MaxDeltaTimeFactor, rSettings.MinNumOfIterations, rSettings.MaxNumOfIterations}; } From f14881b4c49623609c8d1218b24eeaed7e69fe7c Mon Sep 17 00:00:00 2001 From: markelov Date: Thu, 6 Mar 2025 23:35:02 +0100 Subject: [PATCH 15/16] response to Richard's comments --- .../custom_workflows/adaptive_time_incrementor.cpp | 2 +- .../tests/cpp_tests/custom_workflows/test_time_incrementor.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index 7d1715952888..d5a38ed9b892 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -77,7 +77,7 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes { constexpr auto default_min_delta_time = 1.0e-10; - if (rResultantState.Converged()) // it is converged also at the beginning of the cycles + if (rResultantState.Converged()) { if (rResultantState.num_of_iterations < mMinNumOfIterations) { mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); 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 4d5bfe448ec3..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 @@ -10,6 +10,7 @@ // Main authors: Wijtze Pieter Kikstra // Anne van de Graaf // Richard Faasse +// Gennady Markelov // #include "custom_workflows/adaptive_time_incrementor.h" From 3923798c42d0a5581d76a43555305deff2f1ec9e Mon Sep 17 00:00:00 2001 From: markelov Date: Fri, 7 Mar 2025 17:09:06 +0100 Subject: [PATCH 16/16] response to WPK comments 2.0 --- .../adaptive_time_incrementor.cpp | 22 ++++++++++--------- .../adaptive_time_incrementor.h | 6 +++-- .../python_scripts/geomechanics_analysis.py | 6 ++++- 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp index d5a38ed9b892..f410a810e075 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.cpp @@ -25,17 +25,19 @@ AdaptiveTimeIncrementor::AdaptiveTimeIncrementor(double StartTime std::size_t MaxNumOfCycles, double ReductionFactor, double IncreaseFactor, - std::optional MaybeMinDeltaTime, + 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), - mMaybeMinDeltaTime(std::move(MaybeMinDeltaTime)), + mUserMinDeltaTime(std::move(mUserMinDeltaTime)), mMaxDeltaTime(MaxTimeStepFactor * mDeltaTime), mMinNumOfIterations(MinNumOfIterations), mMaxNumOfIterations(MaxNumOfIterations) @@ -75,10 +77,10 @@ double AdaptiveTimeIncrementor::GetIncrement() const { return mDeltaTime; } void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rResultantState) { - constexpr auto default_min_delta_time = 1.0e-10; + 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.Converged()) { if (rResultantState.num_of_iterations < mMinNumOfIterations) { mDeltaTime = std::min(mDeltaTime * mIncreaseFactor, mMaxDeltaTime); } else if (rResultantState.num_of_iterations == mMaxNumOfIterations) { @@ -86,7 +88,7 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes } if (rResultantState.time < mEndTime && mEndTime - (rResultantState.time + mDeltaTime) < - mMaybeMinDeltaTime.value_or(default_min_delta_time)) { + mUserMinDeltaTime.value_or(default_min_delta_time)) { // Up-scaling to reach end_time without small increment mDeltaTime = mEndTime - rResultantState.time; } @@ -97,10 +99,10 @@ void AdaptiveTimeIncrementor::PostTimeStepExecution(const TimeStepEndState& rRes mDeltaTime *= mReductionFactor; } - KRATOS_ERROR_IF(mDeltaTime < mMaybeMinDeltaTime.value_or(default_min_delta_time)) + KRATOS_ERROR_IF(mDeltaTime < mUserMinDeltaTime.value_or(default_min_delta_time)) << "Delta time (" << mDeltaTime << ") is smaller than " - << (mMaybeMinDeltaTime ? "given" : "default") << " minimum allowable value " - << mMaybeMinDeltaTime.value_or(default_min_delta_time) << std::endl; + << (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 911a34ebd991..b1ba9acd8390 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h +++ b/applications/GeoMechanicsApplication/custom_workflows/adaptive_time_incrementor.h @@ -29,7 +29,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public Tim std::size_t MaxNumOfCycles = 10, double ReductionFactor = 0.5, double IncreaseFactor = 2.0, - std::optional MaybeMinDeltaTime = std::nullopt, + std::optional mUserMinDeltaTime = std::nullopt, double MaxTimeStepFactor = 1000.0, std::size_t MinNumOfIterations = 3, std::size_t MaxNumOfIterations = 15); @@ -41,11 +41,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) AdaptiveTimeIncrementor : public Tim private: double mEndTime; + double mTimeSpan; double mDeltaTime; + double mInitialDeltaTime; std::size_t mMaxNumOfCycles; double mReductionFactor; double mIncreaseFactor; - std::optional mMaybeMinDeltaTime; + std::optional mUserMinDeltaTime; double mMaxDeltaTime; std::size_t mMinNumOfIterations; std::size_t mMaxNumOfIterations; diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py index 2778b822539f..fa7126c861c5 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py @@ -34,6 +34,7 @@ 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() @@ -110,6 +111,7 @@ 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) + print("1: ", self.delta_time, t, new_time) self._CheckDeltaTimeSize() # start the new step @@ -156,6 +158,7 @@ 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) @@ -212,7 +215,8 @@ def _GetSimulationName(self): return "GeoMechanics Analysis" def _GetMinDeltaTimeValueOrDefault(self): - return self.min_delta_time if self.min_delta_time is not None else 1.0e-10 + 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