porepy.numerics.time_step_control module
This module contains a routine for iteration-based time-stepping control.
The algorithm is heavily inspired by [1], which was later used in [2].
- [1] Simunek, J., Van Genuchten, M. T., & Sejna, M. (2005). The HYDRUS-1D software package for
simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. University of California-Riverside Research Reports, 3, 1-240.
- [2] Varela, J., Gasda, S. E., Keilegavlen, E., & Nordbotten, J. M. (2021). A
Finite-Volume-Based Module for Unsaturated Poroelasticity. Advanced Modeling with the MATLAB Reservoir Simulation Toolbox.
Algorithm Overview:
Provided recompute_solution = False, the algorithm will adapt the time step based on iterations. If iterations is less than the lower endpoint of the optimal iteration range, then it will increase the time step by a factor iter_relax_factors[1]. If iterations is greater than the upper endpoint of the optimal iteration range it will decrease the time step by a factor iter_relax_factors[0]. Otherwise, iterations lies in the optimal iteration range, and time step remains unchanged.
If recompute_solution = True, then the time step will be reduced by a factor recomp_factor with the hope of achieving convergence in the next time level. The algorithm will keep decreasing the time step unless: (1) the time step is equal to the minimum admissible time step or (2) the number of recomputing attempts has been exhausted. In both cases, an error will be raised.
Now that the algorithm has determined a new time step, it has to ensure three more conditions, (1) the calculated time step cannot be smaller than dt_min, (2) the calculated time step cannot be larger than dt_max, and (3) the time step cannot be too large such that the next time will exceed a scheduled time. These three conditions are implemented in this order of precedence and will override any of the previous calculated time steps.
Algorithm Workflow in Pseudocode:
- INPUT
time_manager // time step control object properly initialized iterations // number of non-linear iterations recompute_solution // boolean flag
- IF time > final simulation time THEN
RETURN None
ENDIF
- IF constant_dt is True THEN
RETURN dt_init
ENDIF
- IF recompute_solution is False THEN
RESET counter that keeps track of number of recomputing attempts IF iterations < lower endpoint of optimal iteration range THEN
DECREASE dt // multiply by an over relaxation factor > 1
- IFELSE iterations > upper endpoint of optimal iteration range THEN
INCREASE dt // multiply by an under relaxation factor < 1
- ELSE
PASS // dt remains unchanged
ENDIF
- ELSE
- IF number of recomputing attempts has not been exhausted THEN
- IF dt is equal to dt_min THEN
RAISE Error // since recomputation will not have any effect
ENDIF SUBTRACT dt from current time // we have to “go back in time” DECREASE dt // multiply by recomputation factor < 1 INCREASE counter that keeps track of number of recomputing attempts
- ELSE
RAISE Error // maximum number of recomputing attempts has been exhausted
ENDIF
ENDIF
- IF dt < dt_min THEN
SET dt = dt_min
ENDIF
- IF dt > dt_max THEN
SET dt = dt_max
ENDIF
- IF time + dt > a scheduled time THEN
SET dt = scheduled time - time
ENDIF
RETURN dt
- class TimeManager(schedule, dt_init, constant_dt=False, dt_min_max=None, iter_max=15, iter_optimal_range=(4, 7), iter_relax_factors=(0.7, 1.3), recomp_factor=0.5, recomp_max=10, print_info=False)[source]
Bases:
objectParent class for iteration-based time-stepping control routine.
- Parameters
schedule (ArrayLike) –
Array-like object containing the target times for the simulation. Unless a constant time step is prescribed, the time-stepping algorithm will adapt the time step so that the scheduled times are guaranteed to be hit.
The schedule must contain minimally two elements, corresponding to the initial and final simulation times. Schedules of size > 2 must contain strictly increasing times. Examples of VALID inputs are: [0, 1], np.array([0, 10, 30, 50]), and [0, 1*pp.HOUR, 3*pp.HOUR]. Examples of INVALID inputs are: [1], [1, 0], and np.array([0, 1, 1, 2]).
If a constant time step is used (constant_dt = True), then the time step (dt_init) is required to be compatible with the scheduled times in schedule. Otherwise, an error will be raised. Examples of VALID inputs for constant_dt = True and dt_init = 2 are: [0, 2] and np.array([0, 4, 6, 10]). Examples of INVALID inputs for constant_dt = True and dt_init = 2 are [0, 3] and np.array([0, 4, 5, 10]).
dt_init (Union[int, float]) – Initial time step. The initial time step is required to be positive and less or equal than the final simulation time.
constant_dt (bool) – Whether to treat the time step as constant or not. If True, then the time-stepping control algorithm is effectively bypassed. The algorithm will NOT adapt the time step in any situation, even if the user attempts to recompute the solution. Nevertheless, the attributes such as scheduled times will still be accessible, provided dt_init and schedule are compatible.
dt_min_max (Optional[tuple[Union[int, float], Union[int, float]]]) –
Minimum and maximum permissible time steps. If None, then the minimum time step is set to 0.1% of the final simulation time and the maximum time step is set to 10% of the final simulation time. If given, then the first and second elements of the tuple corresponds to the minimum and maximum time steps, respectively.
To avoid oscillations and ensure a stable time step adaptation in combination with the relaxation factors, we further require:
dt_min_max[0] * iter_relax_factors[1] < dt_min_max[1], and dt_min_max[1] * iter_relax_factors[0] > dt_min_max[0].
Note that in practical applications, these conditions are usually met.
iter_max (int) – Maximum number of iterations.
iter_optimal_range (tuple[int, int]) – Optimal iteration range. The first and second elements of the tuple correspond to the lower and upper endpoints of the optimal iteration range.
iter_relax_factors (tuple[float, float]) – Relaxation factors. The first and second elements of the tuple corresponds to the under-relaxation and over-relaxation factors, respectively. We require the under-relaxation factor to be strictly lower than one, whereas the over-relaxation factor is required to be strictly greater than one. Note that the values of iter_relax_factors must be selected such that they avoid oscillations in combination with dt_min_max. Refer to the documentation of dt_min_max for the explicit condition to be satisfied.
recomp_factor (float) – Failed-to-converge solution recomputation factor. Factor by which the time step will be multiplied in case the solution must be recomputed. We require recomp_factor to be strictly less than one.
recomp_max (int) – Failed-to-converge maximum recomputation attempts. The maximum allowable number of consecutive recomputation attempts.
not. (print_info. Whether to print on-the-fly time-stepping information or) –
print_info (bool) –
Example
# The following is an example on how to initialize a time-stepping object time_manager = pp.TimeManager(
schedule=[0, 10], dt_init=0.5, dt_min_max=(0.1, 2), iter_max=10, iter_optimal_range=(3, 8), iter_relax_factors=(0.9, 1.1), recomp_factor=0.1, recomp_max=5, print_info=True
) # To inspect the attributes of the object print(time_manager)
- dt_min_max
Min and max time steps.
- Type
tuple[Union[int, float], Union[int, float]]
- iter_optimal_range
Optimal iteration range.
- Type
tuple[int, int]
- iter_relax_factors
Relaxation factors.
- Type
tuple[float, float]
- schedule
List of scheduled times including initial and final times.
- Type
ArrayLike
- compute_time_step(iterations=None, recompute_solution=False)[source]
Determine next time step based on the previous number of iterations.
See also Algorithm Overview and Algorithm Workflow from the module documentation.
- Parameters
iterations (Optional[int]) – Number of non-linear iterations. In time-dependent simulations, this typically represents the number of iterations for a given time step. A warning is raised if iterations is given when recompute_solution = True or constant_dt = True.
recompute_solution (bool) – Whether the solution needs to be recomputed or not. If True, then the time step is multiplied by recomp_factor. If False, then the time step will be tuned accordingly.
- Returns
Next time step if time < final_time. None, otherwise.
- Return type
- static is_schedule_in_simulated_times(schedule, sim_times, rtol=1e-05, atol=1e-08)[source]
Checks if
scheduleis a proper subset ofsim_timesfor given tolerancesReference: https://github.com/numpy/numpy/issues/7784#issuecomment-848036186
- Parameters
- Returns
True if all times in
scheduleintersect the elements ofsim_timeswith relative tolerancertoland absolute tolerance ``atol`. False otherwise.- Return type