A Parallel Scalable Domain Decomposition Preconditioner for Elastic Crack Simulation Using XFEM

In this study, we consider an elastostatic crack problem. As illustrated in Figure 1, the domain $\Omega\subset\mathbb{R}^{2}$ contains an edge crack $\Gamma_{\emph{c}}$, where the internal pressure along the crack surface is $p$. A prescribed displacement $\textbf{\emph{u}}_{0}$ is imposed on part of the boundary $\Gamma_{\emph{u}}$, and the rest of the boundary $\Gamma_{g}=\partial\Omega\backslash\Gamma_{u}$ is subject to a surface force g. The unit normal to the boundaries is denoted as n. The crack surface is distinguished by $\Gamma_{\emph{c}}^{+}$ and $\Gamma_{\emph{c}}^{-}$, such that $\Gamma_{\emph{c}}=\Gamma_{\emph{c}}^{+}\cup\Gamma_{\emph{c}}^{-}$ and with the unit normal to $\Gamma_{\emph{c}}^{+}$ and $\Gamma_{\emph{c}}^{-}$ by $\textbf{\emph{n}}^{+}$ and $\textbf{\emph{n}}^{-}$, respectively. We define $\textbf{\emph{n}}^{-}=-\textbf{\emph{n}}^{+}=\textbf{\emph{n}}$, and $\textit{{p}}^{+}=-\textit{{p}}^{-}=\textbf{{p}}$, where $\textbf{{p}}^{+}=-p\textit{{n}}^{+}$, $\textbf{{p}}^{-}=-p\textbf{{n}}^{-}$. The equilibrium equations and boundary conditions are given by ²⁶

where f is a given body force, $\boldsymbol{\sigma}$ is the stress tensor corresponding to the displacement field u, and $\nabla$ is the gradient operator $\nabla=\left(\partial/\partial x,\partial/\partial y\right)$. For the linear elasticity problem, the relationship between the stress $\boldsymbol{\sigma}$ and strain $\boldsymbol{\varepsilon}$ is

where "trace" is the trace of a tensor, and I is the identity matrix. The material Lamé constants $\mu$ and $\lambda$ are functions of the Young’s modulus E and Poisson’s ratio v,

The solution of (1) lies in the space of admissible displacement field

where $\emph{{H}}^{1}(*)$ is the usual Sobolev space related to the regularity of the solution. A detailed description of the problem in a nonsmooth domain can be found in the work of Grisvard ²⁷. Similarly, the test function space can be defined as

The variational formulation of (1) is to find $\textbf{\emph{u}}\in\emph{{U}}$ such that

where the bilinear form a($\cdot$,$\cdot$) and linear form l($\cdot$) are defined as

To discretize (6), we apply the corrected XFEM²⁸, which uses a level set function to track the position of cracks, the branch enrichment functions to describe the asymptotic field around the crack tips, and the ramp function to handle the blending area. In this study, the computational domain was discretized by a structured Cartesian mesh with quadrilateral elements. As illustrated in Figure 2, the computational domain is initially divided into the crack tip enriched domain $\Omega^{\rm T}$, the blending area $\Omega^{\rm B}$, the Heaviside enriched domain $\Omega^{\rm H}$, and the unenriched domain $\Omega^{\rm U}$. The Heaviside enriched domain $\Omega^{\rm H}$ is the area where all elements are cut by the crack line, excluding the elements that contain the crack tip. The crack tip enriched domain $\Omega^{\rm T}$ is the area within a radius $r_{tip}$ of the crack tip, where all the nodes for each element are enriched by the branch enrichment functions. The blending domain $\Omega^{\rm B}$ is the area around $\Omega^{\rm T}$, where only a part of the nodes for each element are enriched by the branch enrichment functions. The unenriched domain $\Omega^{\rm U}$ is the area without any enriched nodes, such that $\Omega^{\rm U}=\Omega\backslash(\Omega^{\rm T}\cup\Omega^{\rm B}\cup\Omega^{\rm H})$. Note that the Heaviside enriched domain $\Omega^{H}$ may overlap with $\Omega^{T}$ and $\Omega^{B}$, and the nodes in the overlap area are enriched by multiple enrichment functions.

For elements in $\Omega^{\rm U}$, the displacement field is continuous, so the piecewise linear continuous function from the standard FEM can be adopted as the basis function. In the area around cracks, additional basis functions should be adopted to represent the physics of the problem, except for the standard basis function. For the elements in $\Omega^{\rm H}$, the variables cross cracks are discontinuous, and the additional basis function is the Heaviside function, which represents the discontinuity

where $\textbf{\emph{n}}^{-}_{c}$ is the unit normal to the negative side of the crack surface, and $\textbf{\emph{x}}_{c}$ is a point on the crack surface that has the shortest distance to x. For the elements in $\Omega^{\rm T}$, the analytical displacement fields for the combined mode I and II loading can be found in ²⁹. The analytical solutions are not only discontinuous but also includes singularities. The branch enrichment functions can be selected from the basic terms of the analytical solutions.

where $\left(r,\theta\right)$ is the local polar coordinate corresponding to the crack tip. The approximation of the displacement field for the corrected XFEM can be formulated as follows ^{28, 30}

where $\emph{N}_{S}$, $\emph{N}_{H}$, and $\emph{N}_{tip}$ represent the set of all mesh nodes, the nodal subset for elements cut by cracks, and the nodal subset for elements around the crack tips within an enrichment radius $\emph{r}_{tip}$, respectively; $\emph{N}_{i}$, $\emph{N}_{j}$, and $\emph{N}_{k}$ are the standard FEM shape functions; $\emph{u}_{i}$ is the unknown associated with the node set $\emph{N}_{s}$; $\emph{{a}}_{j}$ and $\emph{{b}}_{\alpha}^{k}$ are the unknowns for the enriched nodal subset $\emph{N}_{H}$ and $\emph{N}_{tip}$; $\mathcal{R}(x)$ represents a ramp function that is employed to overcome the problem of blending elements in $\Omega^{\text{B}}$ ³¹. Based on the approximation in (10), the discrete system for (6) has the following formulation:

where F is the force vector, $\textbf{{d}}=(\textit{{u}},\textit{{a}},\textit{{b}})$ is the unknown, and $\textbf{{K}}=(k_{ij})$, known as the stiffness matrix, is an $n\times n$ symmetrical sparse matrix, whose structure is as follows:

where $\textit{{k}}^{uu}$ is the stiffness matrix related to the standard FEM basis, which has the same nonzero structure as the standard FEM; $\textit{{k}}^{aa}$ and $\textit{{k}}^{bb}$ are the stiffness matrices related to the enriched DOFs, which are denser than $\textit{{k}}^{uu}$; and the other parts are the coupling matrices for regular DOFs and enriched DOFs.

In addition, the linear dependency for the local enrichment in each crack tip should be carefully considered, as described in the work of Fries¹⁵. The detailed numerical method for eliminating the linear dependency can be found in Appendix A. Unlike the standard FEM, the matrix K from the corrected XFEM is highly ill-conditioned, and the linear system is difficult to solve ³². An efficient method for solving the ill-conditioned system is to use the preconditioned Krylov subspace with a well-designed preconditioner. The main focus of this study is to design an efficient and scalable preconditioner to solve the large-scale linear system (11).

Because of the singularities of the asymptotic field around the crack tip, the linear system (11) is difficult to solve using classical iterative methods, such as GMRES and CG, especially for large-scale problems in parallel computing. To accelerate the iterative method, preconditioning techniques are required. One of the most popular preconditioning methods in parallel computing is the additive Schwarz method ³³. Solving the linear system by combining the left preconditioned Krylov subspace method with the additive Schwarz method (11) is equivalent to solving the following preconditioned system:

where $\textbf{{M}}_{*}^{-1}$ is the additive Schwarz preconditioner.

The additive Schwarz method is a type of overlapping domain decomposition method. This can be defined in two ways. The first is the algebraic method, which obtains the subdomain matrices $\textbf{\emph{K}}_{i}(i=1,2,...,N)$ and overlaps algebraically based on the matrix K ³⁴. The second method is the geometrical method, which begins by partitioning the physical domain $\Omega$ into nonoverlapping subdomains $\Omega_{i}^{0}$ (i=1, 2, …, N) and then extending each subdomain $\Omega_{i}^{0}$ into an overlapping subdomain $\Omega_{i}^{\delta}$ by adding $\delta$ layers of elements from the neighboring subdomains. For each subdomain, the submatrix $\textbf{\emph{K}}_{i}$ is obtained as follows:

where N is the number of subdomains, $\textbf{{R}}_{i}^{\delta}$ is a restriction operator from the global domain $\Omega$ to the overlapping subdomain $\Omega_{i}^{\delta}$, which is an $n_{i}\times n$ rectangular matrix, where $n$ is the global matrix size and $n_{i}$ is the local matrix size. Only one component in each row of $\textbf{{R}}_{i}^{\delta}$ is set to one, which is related to the DOF belonging to $\Omega_{i}^{\delta}$, and the other components are set to zero. $(\textbf{{R}}_{i}^{\delta})^{T}$ is the extension operator from local to global, which is usually the transpose of the restriction operator. Depending on the method of treating the overlap, there are three types of additive Schwarz preconditioners: classical additive Schwarz preconditioner, restricted additive Schwarz preconditioner, and harmonic additive Schwarz preconditioner. In this study, we adopt the classical additive Schwarz preconditioner (ASM), which is defined as

where $\textbf{{B}}_{i}^{-1}$ is the approximation of the inverse matrix of $\textbf{{K}}_{i}$, which is usually computed by the incomplete Cholesky factorization (ICC).

The numerical results in Section 4.2 indicate that the algebraic additive Schwarz preconditioner and the standard geometrical additive Schwarz preconditioner do not work well for the crack problems considered in this study. The differences in the properties of the subdomains with and without cracks, which can be quite extensive, must be considered when constructing the additive Schwarz preconditioner. In this paper, a geometrical additive Schwarz preconditioner is considered after carefully examining the properties of the crack tip subdomains.

As illustrated in Figure 3, all subdomains are classified into two types: regular subdomains $\Omega_{i}^{r}\ (i\in[1,N_{reg}])$ and crack tip subdomains $\Omega_{j}^{t}\ (j\in[1,N_{tip}])$. The regular subdomains are partitioned based on the geometry of the computational domain $\Omega$, regardless of the crack position. The number of regular subdomains $\textit{N}_{reg}$ is equal to the number of processors $\textit{N}_{p}$. The number of crack tip subdomains $\textit{N}_{tip}$ is equal to the number of crack tips. Thus, $\bigcup\Omega_{j}^{t}=\Omega^{\rm T}\cup\Omega^{\rm B}$ and $\bigcup\Omega_{i}^{r}=\Omega\backslash(\Omega^{\rm T}\cup\Omega^{\rm B})$. For parallel computing, the submatrices from regular subdomains are allocated to each processor locally, whereas the submatrices from the crack tip subdomains may be distributed over arbitrary processors. For convenience of calculation, the submatrices from the crack tip subdomains are aggregated and stored on specific $\textit{N}_{tip}$ processors to compute the inverse matrix of each subdomain according to (15). To decrease the communication between processors, the $N_{tip}$ processors should include the enrichment area around the crack tips as much as possible.

During each iteration of the preconditioned Krylov subspace iterative method, we need to compute $\textbf{{y}}=\textbf{{M}}_{*}^{-1}\textbf{{x}}$, where x is the residual vector or the vector of unknowns defined on the domain $\Omega$. The ASM preconditioner can be expressed as:

This was computed in parallel. That is, on the $i$th processor core, we compute

then obtain the global vector by adding them together

where $\textbf{{y}}_{i}$ includes two parts in the $N_{tip}$ processors, which include a regular subdomain and an enriched crack tip subdomain. The subvector in each processor core with overlap $\textbf{{x}}_{i}^{\delta}$ is defined as

Let $\textbf{{X}}_{i}=\textbf{{R}}_{i}^{\delta}\textbf{{x}}=\textbf{{x}}_{i}^{\delta}$ and $\textbf{{Y}}_{i}=\textbf{{B}}_{i}^{-1}\textbf{{X}}_{i}$, then, $\textbf{{y}}_{i}=(\textbf{{R}}_{i}^{\delta})^{T}\textbf{{Y}}_{i}$. $\textbf{{Y}}_{i}$ is solved by solving the linear system of the equations $\textbf{{K}}_{i}\textbf{{Y}}_{i}=\textbf{{X}}_{i}$ locally. For the specific $N_{tip}$ processors, we solve $\textbf{{K}}_{i}^{reg}\textbf{{Y}}_{i}^{reg}=\textbf{{X}}_{i}^{reg}$ and $\textbf{{K}}_{i}^{tip}\textbf{{Y}}_{i}^{tip}=\textbf{{X}}_{i}^{tip}$, and then $\textit{{y}}_{i}=(\textit{{R}}_{i,reg}^{\delta})^{T}\textbf{{Y}}_{i}^{reg}+(\textit{{R}}_{i,tip}^{\delta})^{T}\textbf{{Y}}_{i}^{tip}$. Here, we distinguish the restriction and extension operators by subscripts $reg$ and $tip$ in the $N_{tip}$ processors. This means that $(\textit{{R}}_{i,reg}^{\delta})^{T}$ and $(\textit{{R}}_{i,tip}^{\delta})^{T}$ are extension operators for regular subdomains and tip subdomains in the $N_{tip}$ processors. A subscript is not needed for the other $N_{p}-N_{tip}$ processors because only one subdomain problem needs to be solved by each processor. Note that the restriction and extension operators with subscript $tip$ deal with all types of DOFs, whereas operators with subscript $reg$ or without any subscript only deal with regular DOFs and Heaviside DOFs. The above algorithm is summarized Algorithm 1.

The main idea of the algorithm is to handle the crack tip subdomain separately. The difference between this method and the standard domain decomposition method (sDDM) is depicted in Figure 4. In sDDM, all the subdomain problems are solved in parallel on each processor, and the communication between the processors is only from the overlaps. In our method, the crack tip submatrix is obtained and stored on specific $N_{tip}$ processors. To reduce communication, a flag in each processor indicates whether the processor has crack tips inside, and the data transfer mainly occurs among processors that include crack tips. If more than two processors share a crack tip area, the one that holds most of the crack tip occupies the subdomain and the problem is solved locally. During each iteration of the preconditioned Krylov subspace method, these $N_{tip}$ processors need to perform more work than sDDM because the linear system from the crack tip area needs to be solved while the other $N_{p}-N_{tip}$ processors are idle and wait for the $N_{tip}$ processors to finish the operation. This affects the synchronization between processors and hence the parallel efficiency. The numerical results in Section 4.7 demonstrate that this effect is very small, as the most time-consuming step is to solve the linear equations from the regular subdomains. It takes very little time to solve the crack tip problem because the size of the crack tip submatrix is usually much smaller than that of the regular subdomain.

In this section, the focus is on the crack-tip area where the difficulty arises. There are two important roles for the branch enrichment functions (9): localization of where the crack is and providing the crack tip area with more information on the analytical solutions¹⁵. The graphical representation of these functions is provided in Figure 8. It is observed that only the first term $\sqrt{r}\sin(\theta/2)$ includes a strong discontinuity, which localizes the position of the crack inside the element. The other three terms, which employ more physics in the solution, are continuous. As discussed in Section 3, the condition number of the stiffness matrix from the XFEM is much larger than that of the standard FEM, which presents significant challenges in problem solving. In applications, only the first few terms ³⁷ are adopted for simplicity, to decrease the number of DOFs and the condition number of the discrete system. The effect of adopting different terms of enrichment functions on the solver and the accuracy of the numerical results is also investigated. The crack model is the same as that described in Section 3.

In Figure 9, the distribution of eigenvalues for different XFEMs is displayed. The number of DOFs ($N_{dof}$) and the condition number of the linear system are also compared. Here, “enr=1” indicates that the crack tip elements are enriched only by the first term of equation (9); “enr=2” indicates that the first two terms are used for enrichment, and so on. It is observed that the maximum eigenvalues remain almost constant, while the minimum eigenvalues decrease quickly with increasing terms (9). The condition number of a linear system is closely related to the distribution of the eigenvalues. A small minimum eigenvalue results in a large condition number. This observation suggests that the main problem arises from the branch enrichment in the crack tip area.

In Figure 10, the L2-norm and SIF errors for the four XFEMs are explored with respect to mesh refinement, whereby with the same mesh size, the full enrichment functions (enr=4) provided the most accurate results. The more the number of enrichment terms used, the better the results obtained. That is, with more enrichment terms, the basis functions in (9) include more information about the analytical solution around the crack tip, and the numerical results are more credible. Table 1 shows the convergence order for the L2-norm, energy-norm, and SIF errors. The results demonstrate that more enrichment terms yield a better convergence order. The XFEM with four enrichment functions can achieve the optimal convergence rate like the standard FEM.

The above discussion indicates that the full enrichment function (enr=4), apart from increasing the condition number of the linear system, also increases the accuracy of the results significantly, and provides the best convergence order. For a linear system that uses part of the branch enrichment functions for the crack problem, conventional iterative solvers can provide good performance. However, the numerical tests demonstrate that these solvers do not converge well for enr=4, unless a good preconditioner is applied to decrease the condition number.

Note that the stiffness matrix resulting from the XFEM discretization of linear elastic problem is symmetric. After applying the boundary conditions and eliminating the linear dependency, the coefficient matrix for the linear system is symmetric positive definite (SPD) ³⁸. For this type of system, the conjugate gradient (CG) algorithm is a good choice for SPD matrices ³⁹. For each subdomain problem, the submatrix is also symmetric, and the complete Cholesky factorization (CC) or incomplete Cholesky factorization (ICC) is adopted as the subdomain solver. Therefore, all cases in this study are solved by CG with an ASM preconditioner for the enr=4 problem, and the subsolver for each subdomain is CC or ICC with $l$ levels of fill-ins, abbreviated as ICC($l$).

In this section, the performances of the different subdomain strategies introduced in Section 2.2 are compared. The branch crack model shown in Figure 13a is employed for the discussion. Young’s modulus $E=10^{4}$ MPa, and Poisson’s ratio $\nu=0.30$. Initially, the condition number of the stiffness matrix is studied before and after preconditioning with a small mesh, as it is difficult to compute the condition number for large size matrices. The mesh size of 26$\times$26 is partitioned into 5$\times$5 subdomains. The overlapping size is 2. There are 1352 regular DOFs, 1088 enriched DOFs, and the assembled stiffness matrix size is 2440$\times$2440. According to the results in Table 2, it is observed that the condition number of the stiffness matrix is very large without the preconditioner. For the preconditioner with strategy S0, improvements do not impact the magnitude of the condition number. Comparing the results of S1 and S2, the minimum singular value and condition number with strategy S2 are noticeably smaller than those of S1. As the magnitude of the condition number decreases from $10^{11}$ to $10^{5}$, the domain decomposition method S2 significantly decreases the condition number.

To compare the performance of the three subdomain strategies in parallel computing, two different meshes with mesh sizes of $1000\times 1000$ and $5000\times 5000$ are adopted for this discussion. The corresponding number of nodes is 1 million and 25 million, respectively. The number of processors used to solve the two problems is 24 and 1024, respectively. The overlapping size is 4 for S0, and 2 for S1 and S2. The subdomain solver of each processor for S0 and S1 is ICC(9). For S2, the subdomain solver is the complete Cholesky factorization (CC) for the tip subdomains and ICC(9) for regular subdomains. As displayed in Table 3, for the linear elastic crack problem, strategy S0 with the preconditioned CG method converges very slowly. When the mesh is refined to 25 million nodes, more than 4500 iterations are required for convergence. The geometrical additive Schwarz preconditioner with strategies S1 and S2 can reduce the number of iterations significantly. Compared with S1, S2 requires fewer iterations and elapsed time as it handles the crack tip subdomain separately, where the difficulties originate. Therefore, the domain decomposition preconditioner from S2 has better performance than the algebraic domain decomposition method S0 and the traditional geometrical domain decomposition method S1.

A desirable property of an algorithm is to maintain a similar convergence rate under different mesh resolutions ^{40, 41} referred to as mesh-independent convergence. In this section, the mesh-independent convergence of the additive Schwarz preconditioner is discussed by comparing the number of iterations and the norm of errors, including SIF-error, L2-norm, and energy-norm, as introduced in Section 3. The crack model remains the same as that of Section 4.2. The mesh is refined from 512$\times$512 to 4096$\times$4096, such that the mesh resolution changes from $h$ to $h/8$, where $h$ is the element size for the coarsest mesh. The number of processor cores for parallel computing is 192, and the overlapping size ranges from 1 to 8 to maintain the thickness of the overlap region between processors constant. Here, the overlapping size for the regular subdomain and crack tip subdomain is the same and the subsolver was the complete Cholesky factorization for both. To compute the norm of errors, the results of the finer mesh size $8192\times 8192$ were adopted as reference. The $u_{i},\varepsilon_{i},\sigma_{i},K_{I}$ in equations (20,21, and 22) without superscript $h$ are the results of the finer mesh.

The convergence performance is presented in Table 4. Clearly, the number of iterations changes slightly with mesh refinement. The SIF error, L2-norm, and energy norm become increasingly smaller. This demonstrates that the performance of the domain decomposition method in this study is independent of mesh resolution. It should be noted that the subsolver in Table 4 is the complete Cholesky factorization instead of the incomplete Cholesky factorization, as our numerical experiments indicate that the subsolver ICC($l$) cannot achieve mesh-independent convergence. Although the convergence of complete Cholesky factorization is mesh-independent, it is too time-consuming and unsuitable for large-scale parallel computing. The subsolver for the regular subdomains in the following discussions is also ICC($l$).

In standard FEM, the DOFs are usually ordered field by field (FbF), which orders one field DOFs and then another. For the corrected XFEM, the regular DOFs are ordered first and then the crack-tip enriched DOFs and Heaviside DOFs. Figure 11a displays the structure of the stiffness matrix. In parallel computing, DOFs are usually ordered block by block (BbB), such that each processor orders local DOFs sequentially, and the starting index of processor $i$ is numbered after processor $i-1$ regardless of the type of DOFs. Figure 11b displays the structure of the stiffness matrix. It is observed that the bandwidths from the two orders are still too large for subdomain solvers such as ICC factorization ⁴². Additional reordering techniques are required to decrease the bandwidth and improve the performance of the subsolver. The matrix reordering methods include ND, 1WD, RCM, and QMD ^{43, 44}. The structure of the stiffness matrix with these four reordering methods is displayed in Figure 11c$\sim$11f. It is observed that RCM can considerably decrease the bandwidth. The performance of the reordering techniques in parallel computing is investigated next. The crack model remains the same as in Section 4.2 and the preconditioner is the left ASM.

In Table 5, the number of iterations and elapsed time of the solver are listed for the mesh with 1 million nodes and 25 million nodes. The DOFs in each processor are ordered using the BbB method by default. If the reordering method is "Natural", this means that there is no additional reordering method. The results show that compared with the natural order, ND has less elapsed time for both cases, while 1WD has no improvement and requires more iterations and elapsed time. RCM requires the least number of iterations and elapsed time among these reordering methods. QMD offers some but not tangible improvements. As RCM performs the best, it will be used as the default matrix reordering technique in the following discussions.

Generally, for large-scale parallel computing, complete Cholesky factorization for a symmetric sparse matrix is very time- and memory-consuming as the complexity is $O(N^{3})$ for an $N\times N$ matrix ⁴⁵. An incomplete Cholesky factorization with $l$ levels of fill-ins, abbreviated as ICC($l$), is an optional method for obtaining the subdomain preconditioner $B_{i}^{-1}$ in (15) because the complexity is $O(N^{2})$. The effect of the ICC fill-ins level on the linear solver is investigated next by comparing the number of iterations and elapsed time. The branch crack model is employed in this section; the mesh is 5000$\times$5000; and the number of unknowns is 50,015,964. The problem is computed by 1024 processors with a left ASM preconditioner; the overlapping sizes for the regular subdomains and crack tip subdomains are 2 and 6, respectively.

The effects of ICC fill-ins level for tip subdomains and regular subdomains are discussed separately. When considering regular subdomains, the subsolver for the crack tip subdomain is fixed as complete Cholesky factorization. When considering crack tip subdomains, the subsolver for the regular subdomain is fixed as ICC(9). The results are summarized in Table 6 and Figure 12. On the left of Table 6, the ICC fill-ins level for the regular subdomain changes from 0 to 16 indicating that the number of iterations continues to decrease. The elapsed time of the solver decreases at first, and finally, there are some increases. Memory usage increases continuously with the increase of the ICC fill-ins level, increasing to a maximum when the subsolver is the complete Cholesky factorization. In Figure 12a, the number of iterations and the elapsed time drop quickly at first and then gradually when the fill-ins level is smaller than 5. When the fill-ins level is larger than 9, the elapsed time begins to increase as the elapsed time per iteration increases continuously, as shown in Figure 12b. This is the reason why $T_{sol}$ increases for large fill-ins level. Compared with the results of complete Cholesky factorization in the last row of Table 6, although the ICC($l$) requires more iterations to converge, the cost for each step is much less, which results in less total time used. On the right part of Table 6, the situation is different. When the ICC fill-ins level for the crack tip subdomain is less than 7, the solver does not converge. With the increase of the ICC fill-ins level, the solver converges and quickly attains the results of CC. This indicates that the crack tip subdomain solver is much more sensitive to the accuracy of matrix factorization because the crack tip submatrix is singular. The iterations, elapsed time, and memory usage are remarkably close to those of complete Cholesky factorization in the last row. As the crack tip submatrix is denser, there is no significant difference between complete Cholesky factorization and incomplete Cholesky factorization with large fill-ins level in terms of time consumption and memory usage. For these reasons, the subsolver for the tip subdomain is set as CC, and that for the regular subdomain is ICC(9), which achieves the best performance.

In this section, the impact of overlapping size, which plays a critical role in the additive Schwarz method is studied. Overlap in the additive Schwarz method represents the communication between processor cores. Generally, a larger overlapping size can decrease the number of iterations because the local processor can obtain more information from neighboring processors. Meanwhile, a larger overlapping size increases the cost of the subsolver for local problems. Therefore, an appropriate overlapping size should be used to achieve the best performance. The crack model is a branch crack; the mesh size is 5000$\times$5000; and the problem is solved using 1024 processors. The subsolver for regular subdomains is ICC(9), and the crack tip subdomain is CC.

Table 7 shows the effect of overlapping size in terms of iterations and elapsed time. Here, the overlapping size for regular subdomains and tip subdomains is discussed separately. On the left part of Table 7, the overlapping size for tip subdomains is fixed as 6, and that for regular subdomains varies from 1 to 7. The number of iterations keeps decreasing continuously, while the elapsed time drops at first and then begins to increase when the overlapping size is larger than 2, as the elapsed time per iteration is always increasing. It is known that each layer of overlapping elements in the geometrical domain decomposition method includes many DOFs, increasing the communication between processors. Therefore, the overlapping size for the regular subdomain can be set to OLP(2). On the right part of Table 7, the overlapping size for regular subdomains is fixed as 2, and that for crack tip subdomains varies from 1 to 7. With the increasing overlapping size, the iterations always decrease and the elapsed time keeps decreasing as well when the overlapping size is smaller than 6. When the overlapping size is 7, there is little increase in elapsed time. The elapsed time per iteration exhibits little change when the overlapping size increases, which can be explained by the crack tip subdomain being usually much smaller than the regular subdomain, leading to much less DOFs from crack tip overlap. The main communication during iterations mainly comes from regular subdomains. The overlapping size for crack tip subdomain has an effect on the number of iterations but not much on the elapsed time per iteration. As summarized, the overlapping size for regular subdomain is chosen as OLP(2) and for the crack tip subdomain, it is OLP(6).

In this section, the parallel scalability of the proposed algorithm is discussed. As shown in Figure 13, a square domain of size 10 m$\times$10 m with cracks inside is considered. The first crack model is a branch crack with three crack tips. The second crack model has 16 cracks in the domain and 32 crack tips. The problem type is plain stress. Young’s modulus $E=10000$ MPa, and Poisson’s ratio $\nu=0.30$. The boundaries of the domain are fixed with zero displacement. The crack surface is compressed with 1.0 MPa inner water pressure. The relative convergence tolerance was $10^{-6}$. The submatrix in each processor for all cases was reordered by RCM. The numerical results for the two crack models are illustrated in Figure 14, which plots the Von-Misses stress contour. It is observed that the high stress is concentrated in the crack tip area. The results agree with the theoretical analysis that singularities near crack tips generate a large stress around. In the following, parallel scalability is investigated by changing the subdomain solver and the number of processors. Based on the discussions above, the subsolvers ICC(8), ICC(9), and CC are considered for the two cases. The overlapping size for regular subdomains is 2, and crack tip subdomains is 6.

For the branch crack model, the mesh size is $5000\times 5000$. After XFEM discretization, there were $5\times 10^{7}$ regular DOFs and 15964 enriched DOFs. The number of processor cores varied from 512, 1024, 1536, to 2048. The numerical results are summarized in Table 8. The number of iterations and elapsed time decrease with an increase of the ICC fill-ins level. When the subdomain solver is set to CC, the iterations decrease to a minimum, whereas the solver time becomes much longer than that of ICC($l$). As for scalability with respect to the number of processors, the iterations increase continuously, and the elapsed time decreases significantly. In the last column for $T_{tip}/T_{reg}$, the subsolver time for crack tip subdomains occupies a relatively small percentage, and has a small impact on the scalability. The speedup and parallel efficiencies are displayed in Figure 15. The speedup with the subdomain solver ICC($l$) keeps increasing but is smaller than the ideal speedup, and the efficiency keeps decreasing with the increase of processors. The situation is different for the subdomain solver CC. The speedup and efficiency both exceed the ideal speedup and ideal efficiency while the rate of increase continues to decrease. This is because the complexity of ICC(l) has a linear dependence on the number of processors, whereas for subsolver CC, the dependence is super-linear.

A more complex model is considered in Figure 13b to test the strong scalability of the proposed algorithm. The mesh is refined to $10000\times 10000$, which has $2\times 10^{8}$ regular DOFs and 123216 enriched DOFs after XFEM discretization. Based on the discussion so far, this large-scale problem is exceedingly difficult to solve because the condition number of the linear system is large (larger than $10^{20}$). Generally, the difficulty in solving a crack problem increases with mesh refinement and the number of cracks. If the domain has more cracks, the linear system includes more singularities, which makes the stiffness matrix more ill-conditioned. In this case, the subdomain solvers were ICC(8), ICC(9), and CC with different numbers of processors: 2048, 4096, 6144, and 8192. The numerical results are presented in Table 9. The results indicate that the basic principle of the subdomain solver and the scalability are the same as that of the branch crack model mentioned above. The difference is that the 16-cracks model requires more iterations and elapsed time to achieve convergence. In the last column for $T_{tip}/T_{reg}$, the subsolver time for the crack tip subproblem is considerable when the number of processors is 8192, as a result of the reduced parallel efficiency. The speedup and efficiency are displayed in Figure 16. Comparing Figure 15 and Figure 16, the speedup and efficiency for the 16-cracks model are lower than those of the branch crack model with the same number of processors and ICC fill-ins level.

The parallel analysis for the two crack models demonstrates that the proposed parallel domain decomposition preconditioner is scalable in terms of the number of processors. For a model with more cracks and refined mesh, the algorithm retains scalability. Despite the increase in the number of iterations and elapsed time to converge, the parallel efficiency remains at more than 70%. Currently, the crack tip subproblem is solved on the part of processor cores, which decreases the parallel efficiency when the number of processors is large although the crack tip solver time is only a small percentage of that of regular subdomains. Such a restriction will be removed in future studies by solving the crack tip subproblem, employing all processors in parallel.

The authors declare no potential conflict of interests.