Time parallelization scheme with an adaptive time step size for solving stiff initial value problems

: In this paper, we introduce a practical strategy to select an adaptive time step size suitable for the parareal algorithm designed to parallelize a numerical scheme for solving stiff initial value problems. For the adaptive time step size, a technique to detect stiffness of a given system is first considered since the step size will be chosen according to the extent of stiffness. Finally, the stiffness detection technique is applied to an initial prediction step of the parareal algorithm, and select an adaptive step size to each time interval according to the stiffness. Several numerical experiments demonstrate the efficiency of the proposed method.


Introduction
We consider numerical techniques to solve sti initial value problems (IVPs) given by where f has continuously bounded partial derivatives up to required order for the developed numerical method. The sti systems are broadly classi ed into two categories -one is to have sti components in just a given system and the other is to have the components in both the system and its solution. In the rst case, solutions of the system behaves smoothly as time is increasing so it can be easily solved by any implicit scheme with an appropriate step size. On the other hand, the solutions of the given system in the second case have sti components expressed as irregularities or sharp fronts in some or whole time intervals. In the intervals, we carefully handle a numerical scheme since the solutions are very rapidly changed. Most interesting research topics, induced from the real applications such as uid dynamics, molecular dynamics, plasma or other physics, are related with the second case. There are lots of numerical strategies to nd e cient and accurate solutions of the sti systems. In this paper, we focus on the parallelization scheme to nd the e cient solutions of the sti IVPs. Time parallelization scheme has received a lot of attention over the past few years and several parallelization schemes have been proposed [1][2][3]. Especially in 2001, a new algorithm which was named parareal algorithm for the solution of time dependent di erential equations in parallel was introduced [4]. It can be de ned by y k n+ = G(t n+ , t n , y k n ) + F(t n+ , t n , y k− n ) − G(t n+ , t n , y k− n ) where the subscript n refers to the time subdomain number, the superscript k refers to the iteration number. F represents a ne propagator, that is, a more accurate solution on a ne grid in time interval [t n , t n+ ] with an initial value y k− n . G represents a coarse propagator, a less accurate approximation in a coarser grid. Note that the F-propagator determines the overall accuracy of the parareal method, whereas the convergence order of the method is decided by the order of the G-propagator and the number of iterations used when it is coupled with a su ciently accurate F-propagator [5,6]. Unfortunately, the traditional parareal scheme has the low parallel e ciency which is bounded by 1/K, where K is the parareal iterations needed to converge to the desired accuracy. In most case, 2 or more iterations are needed, so the e ciency of the traditional parareal scheme is less than 50 percents and even worst in practice. To hurdle this drawback, several advanced parareal techniques based on the deferred correction (DC) methods have been recently introduced [2,5,6], in which DC strategies are utilized within the parareal iteration for the F-propagator by using one or few DC iterations during each parareal iteration.
In the parareal algorithm, there is an important assumption that there are an in nite number of processors to use, so each processor is assigned to each di erent time interval with a uniform time step size. However, only a few nite number of processors can be provided in practice. Even if an in nite number of processors are provided, it is not e cient to assign the uniform step size on each processor without any consideration on the property of the problem, especially for sti systems or partial di erential equations (PDEs) having sharp front. That is, it is more e cient that a larger time step size is assigned for smooth or non-sti regions in solutions, while a smaller time step size is needed for shock or sti regions. Therefore, the usage of adaptive step sizes is very important issue to improve the e ciency of the parareal algorithm. Related to this issue, many researchers have attempted to nd a suitable way to automatically detect sti ness [7].
The aims of this paper are to introduce a criterion to detect sti ness and to develop a scheme for nding an adaptive time step size according to the extent of sti ness in each interval. First of all, for the given system, we need to split the sti and non-sti parts in a given time domain. There are various ways to detect sti ness. For simplicity's sake, we examine a gradient ratio of a given system to split the sti and non-sti parts. Once the sti regions is detected, the corresponding step size should be automatically controlled. So, the time intervals in sti regions should be gradually shrank depending on the extent of the sti ness, while those in non-sti regions are comparatively stretched. Based on these processes, an appropriate time step size for the parareal algorithm is chosen in the sequential step depending on the sti ness at each time interval. Note that in the traditional parareal algorithm, a G-propagator approximates initial values for all time intervals sequentially, with having a uniform step size at the initialization step.
Additionally, a theoretical analysis of the parareal algorithm shows that the stability of the method depends on the choice of G-propagator [5,8]. Especially, for solving highly sti problem, G-propagator should be satis ed an L-stability. Also, each time interval is determined in the initialization step with a G-propagator, the computational cost for G-propagator should be small enough. Overall, Backward Euler (BE) method will be a good candidate for the G-propagator, since BE is unconditionally stable, its computational costs is relatively small and it has L-stability [9], where it can unconditionally ful ll the stability condition of the parareal methods with less computational costs. This paper is organized as follows. In Sec. 2, we brie y describe the original parareal technique and the improved parareal algorithms based on the original one. In Sec. 3, we introduce several parameters for detecting degree of sti ness in each interval and discuss a strategy to select adaptive time step size using sti ness detection to improve the overall e ciency of the parareal algorithm. In Sec. 4, preliminary numerical results are presented to show the e ciency of the proposed scheme. Finally in Sec. 5, future research directions are provided.

Parareal method
In this section, we brie y review the original parareal algorithm and the improved parareal algorithms based on the traditional algorithm.

. Parareal algorithm
As in general parareal algorithm, we assume the time interval [ , T] is divided into N p intervals with each interval being assigned to a di erent processor denoted processors P through P N p . On each interval, the parareal method iteratively computes a succession of approximation y k n+ ≈ y(t n+ ), where k denotes the iteration number. It is de ned using two propagation operators G(t n+ , t n , y n ) and F(t n+ , t n , y n ), for which the propagators search a solution from t n to t n+ using an initial value y n . The G(t n+ , t n , y n ) operator (denoted G) provides a rough approximation of y(t n+ ), the solutions of Eq. (1) with given initial conditions, whereas the F(t n+ , t n , y n ) operator (denoted F) typically gives a highly accurate approximation of y(t n+ ) on the ne discretization of time interval [t n , t n+ ]. Note that typically the G propagator is computationally less expensive than the F propagator, that is, the G propagator is usually a lower order method or computed on a much coarser discretization, while the F propagator is a higher-order method on a ner discretization. So, the parareal method is convergent to a solution of the F propagator applied in serial.
The parareal method begins by sequentially computing y n for n = , . . . , N p , using G propagator, Once each processor has a value y n , the processors can in parallel compute the approximation F(t n+ , t n , y n ).
The parareal algorithm then computes the serial correction step for n = , . . . , N p , The method proceeds iteratively alternating between the parallel computation of F(t n+ , t n , y k n ) and the serial computation of Eq. (4).

. Improved parareal methods
In this subsection, we brie y introduce improved versions of parareal methods to develop for overcoming limitations of the original parareal algorithm.
Although the original parareal algorithm enables us to parallelize numerical algorithms for solving initial value problems, its e ciency is controversial since the low parallel e ciency which is bounded by 1/K, where K is the parareal iterations needed to converge to the desired accuracy. Note that K must be at least 2, the e ciency of the original parareal scheme is less than 50 percents and even worse in practice. To improve the low parallel e ciency for the parareal algorithm, several improved algorithms are developed [2,5,11], in which various deferred correction (DC) techniques are embedded into the parareal framework. In [2,11], a hybrid parareal spectral deferred correction method was introduced in which spectral deferred correction (SDC) strategies are utilized within the parareal iteration, as a F-propagator. Also, in [5], two di erent deferred correction schemes, modi ed DC technique combined with Backward Euler (BE) method and Krylov deferred correction (KDC) [10], are used for G and F propagators in the parareal framework, similar to the hybrid parareal spectral deferred correction method. Commonly in [2,5,11], instead of directly using the SDC scheme or KDC scheme requiring several iterations (SDC sweeps in SDC or Newton-Krylov iterations in KDC) in serial, the F-propagator in each parareal iteration performs one or a few SDC sweeps or Newton-Krylov iterations on the solution from the previous parareal iteration. As the parareal iterations converge, the F solution still converges to the high-accuracy SDC or KDC solution.
The advantage of these schemes is that the F-propagator becomes much cheaper by combining the parareal iterations and DC iterations, compared to a full accurate solver. Because the DC iterations (SDC sweeps in SDC or Newton-Krylov iterations in KDC) are overlapped with the parareal iteration, so the hybrid parareal schemes can unite the two di erent iterations (DC and parareal iterations) [2, p. 281]. Typically, the original e ciency is bounded by K, where K is the number of iterations for the parallel iterations to converge. However, the parallel e ciency of the hybrid parareal SDC or KDC is about K s K, where K s is the number of iterations required of the serial SDC or KDC method to converge to a given tolerance. Note that the hybrid parareal SDC method is a reasonably good choice for non-sti systems and parareal KDC method is suitable for sti systems since KDC was developed to overcome the limitation of SDC for sti ness.

Algorithm
In the parareal algorithm, all processors are typically initialized by using the coarse propagator in a serial way to yield a low accuracy initial condition on each interval which is assigned to each processor. So, all processors except the rst one are idle until passed an initial condition from the previous processor and the idle time is inevasible.

. Sti ness detection
We now describe a parameter which dictates whether there exists a sti ness of a given system in a given time interval. Since sti ness implies that the given system has two di erent time scales. That is, a rate of maximum and minimum value of derivatives for the system in some interval is quite big, then we prescribe that it is sti in the interval. In this respect, this can be easily done by introducing the following measure related to the gradient for the given problem : where n is the dimension of the given system. Note that we restrict the sti ness to the case when sti ness is on the problem and solution. That is, we need to check the change rate of κ(t) since it says whether any big change exists between the previous interval and the current one. This allows us to approximate a conditioning parameter γ(h) in each interval as follows: where h = t m+ − t m . We can easily check that γ(h) goes to when there is a big di erence between κ(t m ) and κ(t m+ ) due to big change in the interval [t m , t m+ ]. Also, κ(t m ) and κ(t m+ ) have quite similar values, γ(h) goes to . That is, the parameter γ(h) can be used to detect the sti ness in the interval [t m , t m+ ], so we say it "sti ness ratio". Therefore, the sti ness ratio γ(h) is used to determine a time step size h according to the degree of the sti ness, since the ratio γ(h) is relatively small when solutions in a given time interval are rapidly changed and γ(h) is large when a change rate of solutions is increasingly small. Note that "the change rate of solutions is small" means that the solutions are smooth in that interval, so the time step size is allowed to be large. Hence the time step size can be selected adaptively small and increasingly large when the sti ness ratio γ(h) is small and relatively large, respectively.

. Adaptive step size selection
For a code implementation, sti ness criteria to choose time step sizes are needed to be set up. When the ratio γ(h) is close to , the step size is expanded and when the ratio goes to , the step size should be shortened. For the sake of simplicity, we amplify the step size twice when the sti ness ration is , which means that there is no di erence of κ(t) in an interval. In addition, the step size is reduced by half when the sti ness ratio is halved. Using these conditions, we simply set up the following criterion to choose a new time step size h new as follows: where h max and h min are constant factors to avoid too fast increase and decrease of the time step, respectively. Based on the parameters and the discussion above, we get the following algorithm to select new step size as follows: . Parareal algorithm with the proposed adaptive step size controller Based on the new step size controller with the sti ness detection discussed in the previous subsections, we present the enhanced parareal methods with adaptive step size. First of all, we assume the time interval of interest [ , T] is divided into N uniform intervals, and each interval [t i , t i+ ] is assigned to a corresponding processor P i . Note that y k i denotes the approximation after the k-th parareal iteration at the i-th node t i .

Predictor Step Decide initial values and step sizes in a serial way
Starting with the initial value y and h , -Get the initial value y i and h i from the previous interval [t i− , t i ] -Using G-propagator, calculate the initial approximation y i+ for t = t i+ on Processor P i in serial.
-Using the step size selection technique discussed above, calculate a new step size h new based on κ and γ. -Send the y i+ and h new to the next interval.

Corrector
Step Parallel Iteration (k + step) for k = , . . . , N − -Using a higher order method, compute F(t n+ , t n , y k n ) on the ne grid in parallel. -After the approximation value y k+ i at t i on each processor P i− is calculated, it is sending to the following processor P i as a new initial value for t i+ .
-Using G-propagator with a new initial value y k+ i , calculate the initial approximation y k+ i+ for t = t i+ on Processor P i . -Update y k+ n+ = G(t n+ , t n , y k+ n ) + F(t n+ , t n , y k n ) − G(t n+ , t n , y k n ), where G(t n+ , t n , y k n ) is the approximation from G-propagator.

Numerical results
In this section, preliminary numerical results are presented to examine the convergence behavior and eciency of the enhanced parareal scheme, compared to the standard implementation of the original parareal scheme.

. Robertson example
As the rst example, we solve a classical problem due to Robertson which describes the kinetics of an autocatalytic reaction given by Robertson. Its system consists of a sti system of 3 nonlinear ODEs given by The initial vector y is given by [ , , ] T . In this experiment, we march from t = to t f = with initial step size h = e − . To investigate the e ectiveness of the proposed scheme, we simply experiment with the adaptive mesh selection using Backward Euler (BE) method and compare the accuracy of the proposed scheme that of the existing method (built-in Matlab function -ode15s). Note that in this experiment, BE is used as a test method since BE is employed for G-propagator of the parareal algorithm in this paper. Fig. 1(a) shows that the proposed step size controller can solve the problem and its result is quite close to that from the existing method ode15s. Note that we just plot the second component of the solution set since the sti ness is on the second component. Also the sti ness is near t = , so we expect the step size is chosen relatively small in this region. For this, we plot the time step size by the proposed scheme over the time domain and compare it with that obtained by the existing method. Fig. 1(b) shows that only 21 time steps are needed to reach the nal time with the proposed scheme, while time steps are required with uniform grid used in the original parareal scheme and 30 time steps with the existing scheme using adaptive time step size. Also, as seen in the gure, larger time step sizes are allowed by the proposed technique in non-sti parts as desired.
Now we apply the adaptive mesh selection strategy to the parareal algorithm with KDC as a F-propagator and BE as a G-propagator. Also, for the experiment, the KDC methods with 4 Radau II nodes are employed, and each parareal iteration performs the 2 outer Newton iterations for desired e ciency and the other conditions such as the tolerance for the Newton-Krylov methods or nonlinear solvers are xed for all simulations. We use a reference solution obtained from KDC scheme with Radau II node and full outer Newton iterations.
To examine the convergence behavior of sti parts, only 10 processors are used and corresponding nal time point is .
. In Fig. 2, we plot the error at the nal time (t = 0.00125) versus the parareal iterations with an initial step size . It can be seen that after a certain number of parareal iterations, the error levels reach a certain tolerance level even for sti parts. It also shows that to reach the nal time .
using uniform grid with a step size h = , it requires 10 processors which is the same number of processors needed for the adaptive mesh parareal scheme. Hence, it must be noted that even using adaptive step sizes, the step sizes in sti parts should become small enough.

. Van der Pol problem
This problem, that models the behavior in an electronic circuit, can be described as a system of two equations given by where y(t) = [y (t), y (t)] T ∈ R and f is de ned by f t, y(t) = y (t), ( − y (t))y (t) − y (t) .
with initial condition y = [ , ] T . From the second component of f , it can be seen that the smaller is, the stronger sti ness of problem is. For the test, we take = and initial time step size h = − . To examine the e ectiveness of the proposed scheme, we plot the solution y over the time domain and its corresponding step size h in Fig. 5. The gure shows that the step size is adaptively chosen to be small in sti parts and increasingly large in non-sti parts. The gure also shows that only 591 time steps are needed to reach the nal time with the proposed scheme. Note that regardless of sti ness, the original parareal scheme have to use the uniform step size using appropriate step size suitable for sti components. For example, if the initial step size h = − , then × time steps are required. It is directly related to the computational time and the number of processes needed in parareal scheme. Now we apply the adaptive mesh selection strategy to the parareal algorithm with KDC as a F-propagator and BE as a G-propagator. Also, for the experiment, the KDC methods with several Radau II nodes are employed, and each parareal iteration performs the 2 outer Newton iterations for desired e ciency and the other conditions such as the tolerance for the Newton-Krylov methods or nonlinear solvers are xed for all simulations.
Since only a few processors (less than 100) can be available in current status, we test parareal algorithm with adaptive step sizes only for non-sti parts.
Using 72 processors, we march t = to t f = . with adaptive time step size and plot the adaptive step sizes over the time domain in Fig. 5. The adaptive step size is almost same as seen in Fig. 5.
With the adaptive step size, we generate numerical results from the parareal algorithm with , and Radau II nodes for F-propagator (KDC) to examine the convergence behavior. Note that we calculate a For the experiment, we plot the error based on the reference solutions for the parareal iteration in Fig. 5. It can be seen that the accuracy of the algorithm after convergence depends on the number of Radau IIa collocation nodes in the KDC methods. Note that the KDC methods using p Radau IIa nodes is converging with an approximate order of p − [5,12,13].