[Core] Load stepping feature

kratos/python_scripts/analysis_stage.py

@@ -3,6 +3,8 @@
 from KratosMultiphysics.process_factory import KratosProcessFactory
 from KratosMultiphysics.kratos_utilities import IssueDeprecationWarning
 from KratosMultiphysics.model_parameters_factory import KratosModelParametersFactory
+from KratosMultiphysics.step_controller import Factory as StepControllerFactory
+from KratosMultiphysics.step_controller import StepController, DefaultStepController
 
 class AnalysisStage(object):
     """The base class for the AnalysisStage-classes in the applications
@@ -59,11 +61,107 @@ def RunSolutionLoop(self):
         """This function executes the solution loop of the AnalysisStage
         It can be overridden by derived classes
         """
+
+        list_of_step_controllers: 'list[tuple[KratosMultiphysics.IntervalUtility, StepController]]' = []
+
+        default_step_controller_settings = KratosMultiphysics.Parameters("""{
+            "interval": [0, "End"],
+            "settings":{}
+        }""")
+
+        if self.project_parameters.Has("list_of_step_controllers"):
+            for step_controller_settings in self.project_parameters["list_of_step_controllers"].values():
+                step_controller_settings.ValidateAndAssignDefaults(default_step_controller_settings)
+                list_of_step_controllers.append((KratosMultiphysics.IntervalUtility(step_controller_settings), self._CreateStepController(step_controller_settings["settings"])))
+
+        if len(list_of_step_controllers) == 0:
+            list_of_step_controllers.append((KratosMultiphysics.IntervalUtility(KratosMultiphysics.Parameters("""{"interval":[0.0, "End"]}""")), DefaultStepController(KratosMultiphysics.Parameters("""{}"""))))
+
+        computing_mp: KratosMultiphysics.ModelPart = self._GetSolver().GetComputingModelPart()
+
+        is_converged = False
         while self.KeepAdvancingSolutionLoop():
-            self.time = self._AdvanceTime()
+            time_begin = self.time # current step lower bound
+            # current step upper bound.
+            #   - This also clones a new time step
+            #   - This sets the DELTA_TIME (time_end - time_begin)
+            time_end = self._AdvanceTime()
+
+            # find the appropriate step controller for the given time frame
+            step_controller = None
+            for interval_utility, step_controller in reversed(list_of_step_controllers):
+                if interval_utility.IsInInterval(time_end):
+                    break
+
+            if step_controller is None:
+                raise RuntimeError(f"No step controller is found for the time period [{time_begin}, {time_end}].")
+
+            step_controller.Initialize(time_begin, time_end)
+
+            self.time = time_end
+
+            # store the positions of the nodes first, because Predict may move the mesh as well, therefore, we
+            # need to collect node positions before moving
+            # this will add a cost for the existing behavior, hence this is guarded by the following if block to
+            # not to have additional cost of collecting data if no stepping is used.
+            if not isinstance(step_controller, DefaultStepController):
+                ta_position = KratosMultiphysics.TensorAdaptors.NodePositionTensorAdaptor(computing_mp.Nodes, KratosMultiphysics.Configuration.Current)
+                ta_position.CollectData()
+
+            # first try to solve for the final time.
             self.InitializeSolutionStep()
             self._GetSolver().Predict()
             is_converged = self._GetSolver().SolveSolutionStep()
+
+            if not is_converged:
+                KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, f"Did not converge for time = {self.time}.")
+
+            current_step_controller_time = time_begin
+            while not step_controller.IsCompleted(current_step_controller_time, is_converged):
+                if is_converged:
+                    KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, f"Step at time = [{time_begin}, {self.time}] converged.")
+                    # so the sub-step converged.
+
+                    # first finalize the success full step
+                    self.FinalizeSolutionStep()
+                    self.OutputSolutionStep()
+
+                    # make the next step ready for solving
+                    computing_mp.ProcessInfo[KratosMultiphysics.STEP] += 1
+
+                    # here we will correctly set TIME and DELTA_TIME
+                    # and put current values in the previous time step
+                    time_begin = self.time
+                    self.time = step_controller.GetNextStep(self.time, is_converged)
+                    computing_mp.CloneTimeStep(self.time)
+
+                    # now collect the nodes positions, which may have moved
+                    # can be used at a later time to reset the nodal positions.
+                    ta_position.CollectData()
+                else:
+                    KratosMultiphysics.Logger.PrintInfo(self.__class__.__name__, f"Step at time = {self.time} did not converge.")
+                    # sub_step did not converge
+                    # do not advance in step. get a new sub-step
+                    self.time = step_controller.GetNextStep(time_begin, is_converged)
+                    computing_mp.ProcessInfo[KratosMultiphysics.TIME] = self.time
+                    computing_mp.ProcessInfo[KratosMultiphysics.DELTA_TIME] = self.time - time_begin
+
+                    # here we need to reset the coordinates of the mesh, if someone has used the move_mesh_flag = true
+                    # reset the mesh coordinates
+                    ta_position.StoreData()
+
+                    # here we reset the solution step data of each node to the previous time step values
+                    # so the Predict can does a better prediction again with the new step
+                    computing_mp.OverwriteSolutionStepData(1, 0)
+
+                self.InitializeSolutionStep()
+                self._GetSolver().Predict()
+                is_converged = self._GetSolver().SolveSolutionStep()
+
+                current_step_controller_time = self.time
+
+            # we need the last finalize solution step, because once the Solver solves, and converges, and
+            # [t_begin, t_end] is reached, it will no longer go in to the stepping while loop.
             self.FinalizeSolutionStep()
             self.OutputSolutionStep()
 
@@ -394,6 +492,9 @@ def _CheckDeprecatedOutputProcesses(self, list_of_processes):
                 IssueDeprecationWarning("AnalysisStage", msg.format(process.__class__.__name__))
         return deprecated_output_processes
 
+    def _CreateStepController(self, step_controller_parameters: KratosMultiphysics.Parameters) -> StepController:
+        return StepControllerFactory(step_controller_parameters)
+
     def _GetSimulationName(self):
         """Returns the name of the Simulation
         """

kratos/python_scripts/step_controller.py

@@ -0,0 +1,227 @@
+from abc import ABC, abstractmethod
+import KratosMultiphysics as Kratos
+
+class StepController(ABC):
+    """
+    Abstract base class for controlling the progression of steps in a simulation.
+
+    Methods
+    -------
+    Initialize(time_begin: float, time_end: float) -> None
+        Prepare the controller for stepping, given the initial and final times.
+
+    GetNextStep(current_time: float, is_converged: bool) -> float
+        Determine the next step size based on the current time and convergence status of the current step.
+
+    IsCompleted(current_time: float, is_converged: bool) -> bool
+        Check if the stepping process has been completed.
+    """
+    @abstractmethod
+    def Initialize(self, time_begin: float, time_end: float) -> None:
+        """
+        Initializes the step controller with the specified start and end times for current time step from the solver.
+
+        Args:
+            time_begin (float): The starting time of the current time step from the solver.
+            time_end (float): The ending time of the current time step from the solver.
+
+        Returns:
+            None
+        """
+        pass
+
+    @abstractmethod
+    def GetNextStep(self, current_time: float, is_converged: bool) -> float:
+        """
+        Determines and returns the next step based on the current simulation time and convergence status.
+
+        Args:
+            current_time (float): The current time in the simulation.
+            is_converged (bool): Flag indicating whether the previous step has converged.
+
+        Returns:
+            float: The time for the next step.
+        """
+        pass
+
+    @abstractmethod
+    def IsCompleted(self, current_time: float, is_converged: bool) -> bool:
+        """
+        Called to check if the stepping process is completed.
+
+        Args:
+            current_time (float): The current simulation time.
+            is_converged (bool): Indicates whether the solution has converged at the current step.
+
+        Returns:
+            bool: True if the stepping is completed, otherwise false.
+        """
+        pass
+
+class DefaultStepController(StepController):
+    """
+    DefaultStepController is a step controller implementation for Kratos that always returns the final time as the next step,
+    effectively disabling intermediate stepping.
+
+    Args:
+        parameters (Kratos.Parameters): Configuration parameters for the step controller.
+
+    Methods:
+        Initialize(start_time: float, time_end: float) -> None:
+            Initializes the controller with the end time of the simulation.
+
+        GetNextStep(current_time: float, is_converged: bool) -> float:
+            Returns the end time as the next step, regardless of the current time or convergence status.
+
+        IsCompleted(current_time: float, is_converged: bool) -> bool:
+            Always returns True, indicating that stepping is completed.
+    """
+    def __init__(self, parameters: Kratos.Parameters) -> None:
+        default_parameters = Kratos.Parameters("""{
+            "type": "default_step_controller"
+        }""")
+        parameters.ValidateAndAssignDefaults(default_parameters)
+
+    def Initialize(self, _: float, time_end: float) -> None:
+        self.__time_end = time_end
+
+    def GetNextStep(self, _: float, __: bool) -> float:
+        return self.__time_end
+
+    def IsCompleted(self, _: float, __: bool) -> bool:
+        return True
+
+class GeometricStepController(StepController):
+    """
+    GeometricStepController implements an adaptive time-stepping controller based on geometric progression for stepping in simulations.
+
+    This controller adjusts the time increment (`delta_time`) dynamically based on the convergence or divergence of the solution at each step. It increases the time step after a specified number of successful convergences and decreases it after a failed attempt, within user-defined bounds.
+
+    Class Methods:
+        GetDefaultParameters() -> Kratos.Parameters:
+            Returns the default parameters for the controller, including factors for convergence/divergence, initial/min/max time steps, and thresholds for incrementing or terminating the step size.
+
+    Constructor:
+        __init__(settings: Kratos.Parameters):
+            Initializes the controller with user-provided or default settings, validating and storing parameters for each interval.
+
+    Methods:
+        Initialize(time_begin: float, time_end: float) -> None:
+            Prepares the controller for a new stepping process, resetting counters and setting initial parameters.
+
+        GetNextStep(current_time: float, is_converged: bool) -> float:
+            Determines the next time step based on the convergence status of the current time step. Increments or decrements the time step according to the configured factors and thresholds.
+
+        IsCompleted(current_time: float, is_converged: bool) -> bool:
+            Checks if the stepping process is completed, i.e., the solution has converged and the end time is reached.
+
+        __SetSubSteppingParameters() -> None:
+            Selects and sets the appropriate sub-stepping parameters based on the current interval.
+
+    Attributes:
+        __divergence_factor: Factor by which the time step is reduced after a failed attempt.
+        __convergence_factor: Factor by which the time step is increased after successful attempts.
+        __delta_t_init: Initial time step size.
+        __delta_t_min: Minimum allowed time step size.
+        __delta_t_max: Maximum allowed time step size.
+        __max_number_of_sub_steps: Maximum number of sub-steps allowed.
+        __number_of_successful_attempts_for_increment: Number of successful steps required to increase the time step.
+        __number_of_failed_attempts_for_termination: Number of failed steps allowed before termination.
+        __sub_stepping_iteration: Counter for the current sub-step iteration.
+        __success_full_attempts_count: Counter for successful attempts.
+        __failed_attempts_count: Counter for failed attempts.
+        __delta_time: Current time step size.
+        __time_begin: Start time of the stepping process.
+        __time_end: End time of the stepping process.
+
+    Raises:
+        RuntimeError: If the maximum number of sub-steps or failed attempts is exceeded, or if the minimum allowed time step is reached.
+    """
+    @classmethod
+    def GetDefaultParameters(cls) -> Kratos.Parameters:
+        return Kratos.Parameters("""{
+            "type"                                       : "geometric_step_controller",
+            "divergence_factor"                          : 0.25,
+            "convergence_factor"                         : 1.5,
+            "delta_t_init"                               : 1.0,
+            "delta_t_min"                                : 1.0,
+            "delta_t_max"                                : 1.0,
+            "max_number_of_sub_steps"                    : 100,
+            "number_of_successful_attempts_for_increment": 2,
+            "number_of_failed_attempts_for_termination"  : 5
+        }""")
+
+    def __init__(self, settings: Kratos.Parameters) -> None:
+        default_parameters = self.GetDefaultParameters()
+
+        settings.ValidateAndAssignDefaults(default_parameters)
+
+        self.__divergence_factor = settings["divergence_factor"].GetDouble()
+        self.__convergence_factor = settings["convergence_factor"].GetDouble()
+        self.__delta_t_init = settings["delta_t_init"].GetDouble()
+        self.__delta_t_min = settings["delta_t_min"].GetDouble()
+        self.__delta_t_max = settings["delta_t_max"].GetDouble()
+        self.__max_number_of_sub_steps = settings["max_number_of_sub_steps"].GetInt()
+        self.__number_of_successful_attempts_for_increment = settings["number_of_successful_attempts_for_increment"].GetInt()
+        self.__number_of_failed_attempts_for_termination = settings["number_of_failed_attempts_for_termination"].GetInt()
+
+    def Initialize(self, time_begin: float, time_end: float) -> None:
+        self.__time_begin = time_begin
+        self.__time_end = time_end
+
+        self.__sub_stepping_iteration = 0
+        self.__success_full_attempts_count = 0
+        self.__failed_attempts_count = 0
+        self.__delta_time = self.__delta_t_init
+
+    def GetNextStep(self, current_time: float, is_converged: bool) -> float:
+
+        self.__sub_stepping_iteration += 1
+        if self.__sub_stepping_iteration >= self.__max_number_of_sub_steps:
+            raise RuntimeError(f"{self.__class__.__name__}: Reached maximum number of Load steps [ max number of Load steps = {self.__max_number_of_sub_steps} ].")
+
+        if is_converged:
+            # the sub-step is converged. So try the next Step
+            self.__failed_attempts_count = 0
+            self.__success_full_attempts_count += 1
+
+            if self.__success_full_attempts_count > self.__number_of_successful_attempts_for_increment:
+                # now we increment the step
+                self.__delta_time = min(self.__convergence_factor * self.__delta_time, self.__delta_t_max)
+                Kratos.Logger.PrintInfo(self.__class__.__name__, f"Incrementing the delta time to {self.__delta_time}")
+
+            return min(self.__time_end, current_time + self.__delta_time)
+        else:
+            # the sub-step is not converged.
+            self.__failed_attempts_count += 1
+            self.__success_full_attempts_count = 0
+
+            if (self.__failed_attempts_count > self.__number_of_failed_attempts_for_termination):
+                raise RuntimeError(f"Reached maximum number of failed attempts [ {self.__failed_attempts_count} / {self.__number_of_failed_attempts_for_termination} ]")
+            else:
+                self.__delta_time = self.__divergence_factor * self.__delta_time
+                Kratos.Logger.PrintInfo(self.__class__.__name__, f"Decrementing the delta time to {self.__delta_time}")
+                if self.__delta_time <= self.__delta_t_min:
+                    raise RuntimeError(f"Reached the minimum allowed delta time [ delta_time = {self.__delta_time}, allowed minimum delta_time = {self.__delta_t_min} ]")
+
+            new_time = current_time + self.__delta_time
+            Kratos.Logger.PrintInfo(self.__class__.__name__, f"Solving Step {self.__sub_stepping_iteration} for time = {new_time:0.6e}")
+            return new_time
+
+    def IsCompleted(self, current_time: float, is_converged: bool) -> bool:
+        return is_converged and abs(current_time - self.__time_end) <= (self.__time_end - self.__time_begin) * 1e-9
+
+def Factory(parameters: Kratos.Parameters) -> StepController:
+    if not parameters.Has("type"):
+        parameters.AddString("type", "default_step_controller")
+
+    step_controller_map = {
+        "default_step_controller": DefaultStepController,
+        "geometric_step_controller": GeometricStepController
+    }
+
+    controller_type = parameters["type"].GetString()
+    if controller_type not in step_controller_map.keys():
+        raise RuntimeError(f"Unsupported controller type = \"{controller_type}\". Followings are supported:\n\t" + "\n\t".join(step_controller_map.keys()))
+
+    return step_controller_map[controller_type](parameters)

kratos/tests/test_KratosCore.py

@@ -103,6 +103,7 @@
 import test_tetrahedral_mesh_orientation_check
 import test_nd_data
 import test_tensor_adaptors
+import test_step_controller
 
 # Import modules required for sequential orchestrator test
 from test_sequential_orchestrator import EmptyAnalysisStage
@@ -233,6 +234,7 @@ def AssembleTestSuites():
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_tetrahedral_mesh_orientation_check.TestTetrahedralMeshOrientationCheck]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_nd_data.TestNDData]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_tensor_adaptors.TestTensorAdaptors]))
+    smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([test_step_controller.TestStepControllers]))
 
 
     if sympy_available: