Mathematical Model of the Dynamics of A Piecewise-Continuous Dislocation Loop: Stiffness Analysis and Selection of Numerical Methods

An analysis of the mathematical model of the dynamics of a piecewise-continuous dislocation loop is carried out. It is shown that it has pathological stiffness. Due to high stiffness and large dimensionalityof the model it has been proved that the currently used numerical Gear method is inefficient in the implementation of this model (which is confirmed by test experiments). The paper provides an overview of numerical methods suitable for solving stiff ODE systems. As consequence of the review, it has been found that numerical methods ParSODE and PETScare most effective in the implementation of the mathematical model of the dynamics of a piecewise-continuous dislocation loop.


I. INTRODUCTION
The sweeping area of anelementary crystallographic slip is usually much smaller than the size of the deformed crystal, and only a small fraction reaches the surface of the crystal. The time of formation of dislocation clusters and ensembles is much smaller than the duration of the deforming impact. Experimental methods are most effective in the study of the deformation structure of metal at some instant. Currently,to study the formation dynamics of dislocation structures, methods of simulation and mathematical modeling are used [1].
In the early 2000s, a team of the scientific school named after L.E. Popov had developed a mathematical model [2], in which the obstacle field of a dislocation or of a different nature is replaced with a homogeneous medium, rendering the same resistance to the moving dislocation as the original obstacle field [3]. In the model [2] the assumption was used that a dislocation loop after closing has a circular shape and retains it upon the expansion. Despite a rather rough assumption, the results of computational experiments with the use of the model [2] are in good agreement with the experimental data [4] and can be used in the evaluation of some dynamic characteristics of dislocation loop formation.
To account for the form-change in dislocation loops,a mathematical model [5] [7] have shown interesting results, also consistent with the experimental data. However, the used in the calculation Gear method, effective in the implementation of the model [2], proved to be ineffective in the implementation of the model [5]. The paper presents an analysis of the mathematical model [5] on stiffness and the review of numerical methods for its implementation.

II. ANALYSIS OF THE MATHEMATICAL MODEL
The mathematical model [5] is presented in the form of a system of ordinary differential equations (ODE), with the dimensionalitytwo times greater in regards to the quantity of the considered directions of its motion in the modeling the expansion of the dislocation loop. Empirically, it has been shown that to achieve good accuracy the dimensionality of the ODE system [5] must be less than 720.
Today, a common and reliable method for determining the stiffness does not exist. One of the most common definitions of stiffness [8] will be used in the paper: the system is called stiff if all eigenvalues of the system satisfy the following conditions: whereN is the dimensionality of the ODE system, S is the stiffness coefficient [9]. The eigenvalues of the ODE system of the mathematical model [5] are functions of the independent variable (time),the stiffness at the integration interval is not constant. To study changes in the stiffness of the Cauchy problem [5], the value of the stiffness coefficient is calculated in different parts of the integration interval. The results of the study onthe changein the stiffness of the Cauchy problem [5] are shown in Fig. 1 for the first emitted dislocation loop by the dislocation source. The parameter values of the mathematical model [5] specific to copper at room temperature have been used in calculations. The eigenvalues, obtained in allcalculations, have a negative value. The value of the stiffness coefficient is large (Fig. 1 a). All eigenvalues have imaginary parts, and the maximum in modulus real part of the eigenvalue varies from the 9 th to the 12 th order (Fig. 1 b). Thus, in accordance with the noted abovedefinition, the Cauchy problem for the mathematical model [5] is stiff. Moreover, since S ≥ O(10 9 ), then according to [9] the model [5] has a pathological stiffness.

III. REVIEW OF NUMERICAL METHODS
From the above mentioned analysis it follows that for the implementation of the mathematical model [5] it is necessary to use special numerical methods effective for problems with large dimensionalityand pathological stiffness. It is obvious, that the implementation of numerical methods must be performed using parallel or distributed computing systems. The authors have analyzed the numerical methods and solvers most appropriate for problem solving of such a level. Below is a brief description of the most popular ones.

A. VODE
VODE (a Variable-coefficient ODE solver) is a package of subprograms for numerical solution ofODE systems of the first order [10]. The package may be used for both stiff and non-stiff systems. VODE has a variable-coefficient and a fixed leading coefficient form of the backward differentiation formula (BDF). Jacobi matrix is considered as a full or a band with the possibility of multiple useunder certain conditions [11]. The implementation of the VODE method in C/C ++ is called CVODE.

B. PVODE
PVODE is a parallel implementation of the CVODE method with preconditioned Krylov iteration methods [11]. The parallelism in PVODE is achieved by distributing the vectors of ODE solution into user-defined segments, and parallelizinga set of vector cores, respectively [12].

C. WRVODE
WRVODE (waveform relaxation VODE) was introduced as an efficient parallel method for solving large sparsely coupled ODE. Method WRVODE works faster than PVODE in the case when there are frequent delays in the calculation (for synchronization and so forth), but slower when there are no delays [13]. Therefore, WRVODE can be recommended for the implementation in GRID systems.WR method is used with other numerical methods. For example, WR method with uses Euler and fixed steps. The peculiarity of this algorithm is a very good scalability (up to 100 units were used in the calculations). In [13] it has been shown that WR method with a implements Euler is carried out faster than PVODE at delays.

D. PETSc
PETSc is a collection of data structures and subprograms for a scalable (parallel) solution of scientific problems described by PDE and ODE. In PETSc, CPUs is supported through MPI, and GPUs through CUDA or OpenCL [14]. To solve stiff ODE systems inPETScthePVODE method is used [15].

E. ParSODE
The basis of the parallel code ParSODE is the method MIRK [16]. ParSODE in most cases is significantly more effective than VODE and RADAU. The resultant advantage in speed is about 3-5 times. However, VODE and RADAU operate more efficiently when it comes to small ODE systems oflow complexity.The conclusion is that ParSODEiseffective for large ODE, and F. LSODE LSODE (Livermore Solver for Ordinary Differential Equations) is designed for solving stiff and non-stiff ODE systems [17]. Adams methods are used in solution of non-stiff systems (predictor-corrector), BDF methods are used in solution of stiff problems. Problems with unstable high-speed transient processes can cause inefficient or unstable operation of the method. In [18] it is noted that for relatively simple ODE with small dimensionality the PDIRK method is more efficient than the LSODE. The fixed-leading-coefficient version of VODE ran slightly faster than LSODE [19].

G. LSODPK
LSODPK is a package of subprograms on FORTRAN, which uses a merge of LSODE and preconditioned Krylov subspace methods [20]. The corrector iteration composed of Newton iteration or one of the four preconditioned Krylov subspace methods are usedto solve stiff systems in LSODPK. This method is a nonmatrix, which eliminates the explicit Jacobian storage. Krylov methods are ideal in solution of large-scale stiff ODE and provide an advantage in accuracy as compared to many modern solvers, in particular with RADAU [21].
H. RK RK (Runge-Kutta) isoneofthefirstmethodsforsolvingODE. In the case of moderate stiffness, the RK method is faster than VODE, especially in the case when using a graphic processor. However, the solver VODE becomes more efficient with an increase in the stiffness [22].
The RK method is clearly inferior to several other methods (solvers) when it comes to solving systems of large dimensionality and stiffness [23]. For example, EPISODE (old solver of stiff ODE systems in which the Gear method is implemented) copes with such ODE systems more efficiently.
Implicit integration algorithms based on BDF methods are used inmany solvers of stiff ODE systems, but they require costly operations of linear algebra on the Jacobian matrix. Furthermore, a complex logic flow results in divergent directions due to different initial conditions, making implicit algorithms unsuitable for operation on GPUs. A stable explicit scheme such astheRunge-Kutta-Chebyshev (RKC) method is an alternative implicit algorithm for solving problems with a moderate stiffness level. In general, RKC provides a more accurate result than LSODE [24].

IV. CONCLUSION
Thus, it has been show in the paper that the mathematical model of the dynamics of a piecewise-continuous dislocation loop [5] is a problem with an ample quantity of non-linear equations (large dimensionality) and pathological stiffness. The carried out analysis of numerical methods for solving stiff ODE systems has allowed to select ParSODE and PETScfor the implementation of the model [5]. The choice is justified by the possibility of scaling the calculations on CPU or GPU for solving large-scale problems, as well as the possibility to be used in the case of very high stiffness.Moreover, libraries are free. It is planned to apply the selected solvers in the implementation of the model [5] on the supercomputer SKIF Cyberia, located in Tomsk State University [25]. It has 564 2xXeon 6500 processors and in 2007 SKIF Cyberia was number 199 to 500 in Top-list. It is also plannedto carry out a more accurate comparison of solvers PETSc and ParSODE at different parameter values of the model [5]before a comprehensive study on the dynamics of the expansion of closed dislocation loops forming a crystallographic shear zone.