A High-Performance Implementation of a Robust Preconditioner for Heterogeneous Problems

TL;DR

This paper introduces an efficient, scalable implementation of the GenEO preconditioner within the DUNE framework, demonstrating its robustness and high performance on large, complex industrial problems with over 200 million degrees of freedom.

Contribution

It presents a high-performance, scalable implementation of the GenEO preconditioner, showcasing its effectiveness on large-scale heterogeneous problems and providing detailed technical insights.

Findings

Successfully scaled to over 15,000 cores

Maintains effectiveness on complex parameter distributions

Outperforms established methods in intractable scenarios

Abstract

We present an efficient implementation of the highly robust and scalable GenEO preconditioner in the high-performance PDE framework DUNE. The GenEO coarse space is constructed by combining low energy solutions of a local generalised eigenproblem using a partition of unity. In this paper we demonstrate both weak and strong scaling for the GenEO solver on over 15,000 cores by solving an industrially motivated problem with over 200 million degrees of freedom. Further, we show that for highly complex parameter distributions arising in certain real-world applications, established methods become intractable while GenEO remains fully effective. The purpose of this paper is two-fold: to demonstrate the robustness and high parallel efficiency of the solver and to document the technical details that are crucial to the efficiency of the code.

Tables

Table 1. Table 1 : Demonstration of performance of different preconditioners for a fixed problem size of 30,000 DOFs and increasing the number of subdomains: Number of CG iterations (it), coarse space dimension ( dim ( V h ) dimension subscript 𝑉 ℎ \dim(V_{h}) ), an estimate of the condition number κ ​ ( 𝑨 ) 𝜅 𝑨 \kappa(\boldsymbol{A}) .

	1-level		GenEO			BoomerAMG
$N_{\text{core}}$	iter.	$\kappa \left(𝑨\right)$	iter.	$\kappa \left(𝑨\right)$	dim($V_H$)	iter.
4	89	79,735	16	10	78	258
8	97	84,023	15	9	126	258
16	107	98,579	16	10	182	257
32	158	226,871	16	9	526	263

Table 2. Table 2 : Parallel performance of the composites application on Archer . Left: Details of the weak scaling test. The number of elements per core was fixed at 2808 2808 2808 . Right: Strong scaling test demonstrating near optimal strong scaling up to at least 11 , 320 11 320 11,320 cores.

$N_{core}$	DOF	iter.	$\kappa \left(A\right)$	$dim\left(V_H\right)$	Time (sec)
480	$6.4\cdot {10}^6$	156	445	5025	734
960	$1.3\cdot {10}^7$	154	421	7840	806
1920	$2.6\cdot {10}^7$	152	322	18752	800
3840	$5.1\cdot {10}^7$	144	287	29444	772
7680	$1.0\cdot {10}^8$	132	303	50930	764
15360	$2.0\cdot {10}^8$	102	245	94527	845

Equations

$$a(v,w)=⟨f,w⟩,\forall w\in V,$$

$$a(v,w)={\int }_{\Omega }\kappa (x)\nabla v(x)\cdot \nabla w(x)dx={\int }_{\Omega }f(x)w(x)dx,\forall w\in V.$$

$$a(v,w)={\int }_{\Omega }C(x)\epsilon (v):\epsilon (w)dx={\int }_{\Omega }f\cdot wdx+{\int }_{\partial \Omega }(\sigma \cdot n)\cdot vdx,\forall w\in V,$$

$$a(v_h,w_h)=⟨f,w_h⟩,\forall w_h\in V_h.$$

$${\Omega }_j^o:=\{x\in {\Omega }_j:\exists j^{\prime }\ne j\mathrm{such}\mathrm{that}x\in {\Omega }_{j^{\prime }}\}.$$

$$⟨R_{j}g,v⟩=⟨g,R_j^{T}v⟩,$$

$${\mathbf{M}}_{\mathrm{AS},1}^{-1}:=\sum _{j=1}^N{\mathbf{R}}_j^T{\mathbf{A}}_j^{-1}{\mathbf{R}}_j.$$

$${\mathbf{M}}_{\mathrm{AS},2}^{-1}:={\mathbf{R}}_H^T{\mathbf{A}}_H^{-1}{\mathbf{R}}_H+\sum _{j=1}^N{\mathbf{R}}_j^T{\mathbf{A}}_j^{-1}{\mathbf{R}}_j.$$

$${\Xi }_j(v):=\sum _{k\in \mathrm{dof}({\Omega }_j)}{\mu }_{j,k}{v}_k{\varphi }_k{|}_{{\Omega }_j},\text{ for any }v\in V_h({\Omega }_j).$$

$$a_{{\Omega }_j}(p,v)=\lambda a_{{\Omega }_j^o}({\Xi }_j(p),{\Xi }_j(v)),\forall v\in V_h({\Omega }_j).$$

$$V_H:=\mathrm{span}\{R_j^T{\Xi }_j(p_k^j):k=1,\dots ,m_j;j=1,\dots ,N\}.$$

$$m_j:=\min \left\{m:{\lambda }_{m+1}^j>\frac{{\delta }_j}{H_j}\right\},$$

$$\kappa ({\mathbf{M}}_{\mathrm{AS},2}^{-1}\mathbf{A})⩽(1+k_0)\left[2+k_0(2k_0+1)\underset{1⩽j⩽N}{\max }\left(1+\frac{H_j}{{\delta }_j}\right)\right].$$

$${\mathbf{A}}_j{p}_k^j={\lambda }_k^j{\mathbf{X}}_{\mathbf{j}}{\mathbf{A}}_j^o{\mathbf{X}}_j{p}_k^j,$$

$$({\mathbf{A}}_H{)}_{i,j}=({\mathbf{R}}_H{\mathbf{A}}_H{\mathbf{R}}_H^T{)}_{i,j}={\phi }_i{\mathbf{A}}_i{\phi }_j.$$

$$({\mathbf{R}}_H{v}_h{)}_i={\phi }_i\cdot v_h.$$

$${\mathbf{R}}_H^T{v}_H=\sum _i{\phi }_i(v_H{)}_i.$$

Full text

1 Introduction

Computer simulations have become a vital tool in science and engineering. The demand for solving PDEs on ever larger domains and increasing accuracy necessitates the use of high performance computers and the implementation of efficient parallel algorithms. When designing parallel algorithms two issues are crucial: robustness and scalability.

1. i)

Robustness: The parameters involved in the PDE affect the performance of the algorithm to a large extent. A frequent issue is a distribution of parameters with a large contrast and that contain jumps at different length scales. In many cases a large contrast leads to very slow convergence or even stops it completely. 2. ii)

Scalability: The immediate scalability of the finite element method is limited as each degree of freedom is coupled with all others.

One approach to achieve scalability in solving partial differential equations are domain decomposition methods (see e.g. [11, 14]), splitting the given domain into multiple subdomains. The solution of the original problem restricted to each subdomain is computed in parallel and the results are combined to form an approximate solution. This is repeated until convergence is reached. The number of these iterations, however, still depends strongly on the number of subdomains involved as well as coefficient variations. Introducing an additional coarse space that covers all subdomains can restore performance for large numbers of subdomains.

This global space can be tailored to specific problem, as in the generalized finite element method [3]. While these methods are applicable to an entire class of parameter distributions, each of these approaches is based on certain assumptions on the parameters, e.g. the parameters vary strongly only in one direction. The GenEO coarse space chosen in this work is a related approach, originally introduced in [13]. It does not require a-priori knowledge of the parameter distribution and is applicable to a wide range of problems making it suitable as a ’black-box’ solver.

In this paper we focus on two different elliptic problems; the Darcy equation describing incompressible flow in a porous medium and the anisotropic linear elasticity equations. For both equations the case of heterogeneous coefficients is of great interest.

Composite materials, which make up over 50% of recent aircraft constructions, are manufactured from carbon fibres and soft resin layers. The large jump in material properties between the layers makes the simulation of these materials challenging. Commercial solvers such as ABAQUS often rely on direct solvers to deal with these jumps [9]. However, the scalability of direct solvers is limited. We will demonstrate that the GenEO approach converges independentally of the contrast in material properties and the number of subdomains.

The paper is structured as follows. We will sketch the construction of the GenEO preconditioner. Then we will discuss how to efficiently implement the solver in the high-performance finite element framework DUNE [4, 5]. Finally we will provide several numerical experiments demonstrating both the robustness and scalability of the solver, including one large-scale industrially motivated example.

2 Problem formulation and variational setting

Let $V$ be a Hilbert space, $a\,:\,V\times V\rightarrow\mathbb{R}$ a symmetric and coercive bilinear form and $f\in V^{\prime}$. We consider the following abstract variational problem. Find $v\in V$ such that

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the duality pairing.

This variational problem is associated with an elliptic boundary value problem on a domain $\Omega\subset\mathbb{R}^{d},$ $d=2,3$ with Dirichlet boundary $\partial\Omega_{D}$. In particular, we focus on the following two examples.

1. i)

**Darcy problem: ** Given material properties $\kappa\in L^{\infty}(\Omega)$, find $v\in V=\{v\in H^{1}(\Omega)\,:\,v|_{\Omega_{D}}=0\}$ such that

[TABLE] 2. ii)

**Linear Elasticity: ** Given material properties $C$, find $v\in V=\{v\in H^{1}(\Omega)^{d}\,:\,v|_{\Omega_{D}}=0\}$ such that

[TABLE]

where $\varepsilon_{ij}(v)=\frac{1}{2}(\partial_{i}v_{j}+\partial_{j}v_{i})$ is the strain, and $\sigma_{ij}(v)=\sum_{k,l=1}^{d}C_{ijkl}\varepsilon_{kl}$ is the stress.

Consider a discretization of the variational problem (1) using finite elements on a mesh $T_{h}$ of $\Omega$ such that $\overline{\Omega}=\cup_{\tau\in T_{h}}\tau.$ Let $V_{h}\subset V$ be a conforming space of finite element functions. Then the discrete form of (1) is: Find $v_{h}\in V_{h}$ such that

[TABLE]

3 The GenEO Preconditioner

In order to leverage the potential of modern high performance computers, parallelization is crucial. The task must be split into pieces that can be computed independently, and communication between processes must be minimized. However, in the finite element method each degree of freedom is indirectly coupled with all others. In the overlapping additive Schwarz method parallelization is achieved by splitting the computational domain $\Omega$ into multiple overlapping overlaps subdomains and solving a local problem on each subdomain. In an iterative procedure results from each subdomain are added on the overlaps using nearest-neighbor communication and an updated local problem is solved taking into account new information from neighboring subdomains.

We generate the overlapping subdomains by starting from a non-overlapping subdivision $\{\Omega_{j}^{{}^{\prime}}\}_{j=1}^{N}$ of $\Omega$. Each $\Omega_{j}^{{}^{\prime}}$ is the union of mesh elements from $\mathcal{T}_{h}$. An arbitrary number of layers of elements can be added to each $\Omega_{j}^{{}^{\prime}}$ by applying the definition multiple times, resulting in overlapping subdomains $\Omega_{j}$. For each subdomain $1\leqslant j\leqslant N$, the overlapping zone is defined as

[TABLE]

We denote the restriction of the function space $V_{h}$ to $\Omega_{j}$ by $V_{h}(\Omega_{j}):=\{v|_{\Omega_{j}}:v\in V_{h}\}$ for each $1\leqslant j\leqslant N$ and the restiction of $V_{h}(\Omega_{j})$ to functions that inherently fulfill a homogeneous Dirichlet condition on their respective subdomain by $V_{h,0}(\Omega_{j})$.

Definition 1

We define the prolongation operator $R_{j}^{T}:V_{h,0}(\Omega_{j})\rightarrow V_{h}$ for each element $v_{j}\in V_{h,0}(\Omega_{j})$ as $R_{j}^{T}v_{j}|_{\Omega_{j}}=v_{j}$ and $R_{j}^{T}v_{j}|_{\Omega\backslash\Omega_{j}}=0.$ The corresponding restriction operator by $R_{j}$ is defined as

[TABLE]

We denote the matrix form of the restriction operators $R_{j}$ by $\boldsymbol{R}_{j}$ and of the system matrix by $\boldsymbol{A}$. Further, we denote the problems restricted to the subdomains by $\boldsymbol{A}_{j}:=\boldsymbol{R}_{j}\boldsymbol{A}\boldsymbol{R}_{j}^{T}$, for all $j=1,\ldots,N$. Then the additive Schwarz preconditioner is given by

[TABLE]

Due to local exchange of information, the number of iterations required tends to increase strongly with the number of subdomains involved. This can be overcome by additionally solving a suitable global coarse problem. The resulting method is referred to as a two-level additive Schwarz method. For the coarse space $V_{H}$, denote the natural embedding by $R_{H}^{T}:V_{H}\rightarrow V_{h}$ and its adjoint by $R_{H}$.

Definition 2 (Two-level Additive Schwarz)

Denote the problem restricted to the coarse space by $\mathbf{A}_{H}:=\boldsymbol{R}_{H}\boldsymbol{A}\boldsymbol{R}_{H}^{T},$ for all $j=1,\ldots,N$. Then the two-level preconditioner is given by:

[TABLE]

The analysis framework from [14] allows for upper and lower bounds on the condition number of the preconditioned system. The largest eigenvalue of the preconditioned system can be bounded using the maximum number of subdomains that cover each point $k_{0}$. Clearly $k_{0}$ can easily be controlled by constructing a reasonable domain decomposition. The bound on the smallest eigenvalue depends on the stable splitting constant. Thus, this constant should ideally be small. The GenEO coarse space was designed to minmise the stable-splitting constant, it was first introduced in [13], followed by a full theoretical analysis in [12].

A key ingredient in the GenEO coarse space is the partition of unity, which allows ’stitching together’ the local basis results on each subdomain to form a suitable basis of the entire domain.

Definition 3 (Partition of unity)

Given weights $\mu_{j,k}\in[0,1]$ with $\sum_{1\leq j\leq N}\mu_{j,k}=1$, the associated partition of unity operator for each subdomain $1\leqslant j\leqslant N$ is defined by

[TABLE]

Here we denote the value of $v$ at degree of freedom $k$ as $v_{k}$, and by $\phi_{k}$ the basis function associated to the same degree of freedom.

Figure 1 shows two examples of partitions of unity. It should be noted that implementing a smooth partition of unity for arbitrary subdomain arrangements can be challenging and in our tests we have not found this to have a large influence on the efficiency of the preconditioner. Using the partition of unity operators defined above, we can now state the definition of the GenEO coarse space. The space is constructed from the eigenvectors of a specific generalized eigenproblem representing the inequality required by the stable splitting.

For each subdomain $j=1,\ldots,N$, we define the generalized eigenproblem: Find $p\in V_{h}(\Omega_{j})$ such that

[TABLE]

Note that the eigenproblems are local to their respective subdomain $\Omega_{j}$, i.e. they can be computed in parallel. To use them as a global basis they need to be extended to the entire domain using the partition of unity operators.

Definition 4 (GenEO coarse space)

For each subdomain $j=1,\ldots,N$, let $(p_{k}^{j})^{m_{j}}_{k=1}$ be the eigenfunctions from the eigenproblem in (5) corresponding to the $m_{j}$ smallest eigenvalues. Then the GenEO coarse space is defined as

[TABLE]

In [12] the following bound on the condition number of the matrix has been shown.

Theorem 3.1

For all $1\leqslant j\leqslant N$, let the number of eigenvectors chosen in each subdomain be

[TABLE]

where $\delta_{j}$ is a measure of the width of the overlap $\Omega_{j}^{o}$ and $H_{j}=\operatorname{diam}(\Omega_{j})$. Then,

[TABLE]

Thus the algorithm presented is provably independent of the number of subdomains and coefficient variations in the problem. Note that the choice of threshold by which the eigenvectors are selected is not unique. The threshold can be scaled to vary the condition number of the preconditioned system while retaining robustness. This allows us to control the total number of preconditioned CG iterations.

4 HPC Implementation of GenEO in Modern PDE Frameworks

When implementing the GenEO preconditioner in a PDE framework, the primary goal is to preserve the beneficial properties offered by its theoretical construction, namely:

1. i)

**High parallel scalability: ** Since the condition bound in Thm. 3.1 is independent of the number of subdomains we expect the implementation to yield high parallel scalability. The solution of the eigenproblems parallelizes trivially. However, care has to be taken when it comes to the communication necessary to set up the coarse matrix. 2. ii)

**Robustness with respect to problem parameters: ** While this is an inherent property of the preconditioner, some care is required in implementing the Dirichlet boundary conditions. 3. iii)

**Applicability to various types of PDEs: ** The theoretical framework only requires a symmetric positive definite bilinear form as in (1). This flexibility can be preserved in any numerical framework that is based on abstract bilinear forms. This is the case for many modern PDE frameworks, e.g. FEniCS [1], DUNE [4], or deal.ii [2].

In this section, we present a new implementation of the GenEO coarse space and preconditioner within DUNE (Distributed and Unified Numerics Environment), which fulfills these properties. This serves as a reference for the implementation, which is freely available as part of the dune-pdelab module [6] since version 2.6, as well as a general guideline for future implementations in other software packages. DUNE is a generic package that provides the user with key ingredients for solving any FEM problem. As an open source framework written using modern C++ programming techniques, it allows for modularity and reusability while providing HPC grade performance.

4.1 Prerequisites

Many of the components required to implement a two-level Schwarz method already exist within DUNE. In particular, we use the PDELab discretization module’s functionality to assemble stiffness matrices based on bilinear forms and for efficient communication across overlapping subdomains. The GenEO basis functions have support not restricted to individual elements, which makes the existing high-level components of PDELab unsuited for storing the coarse space. As part of this project, components facilitating such coarse spaces were fully integrated within the framework. Further, an efficient sequential solver for generalized eigenproblems is needed. Here, we choose ARPACK [8].

4.2 General Structure

The implementation in PDELab closely follows the structure of the previous section. All mathematical objects are represented as individual classes (see Fig. 2(b)). This separation of concerns leads to an easy to understand and well-structured code. Further, components are easily interchangeable when constructing related methods. In particular, the intricate process of constructing a global coarse space from per-subdomain basis functions is entirely contained in the class SubdomainProjectedCoarseSpace. Thus, the GenEO basis can easily be replaced by a different approach, as only the local basis functions need to be defined on that level.

4.3 Discrete Basis

To calculate GenEO basis functions we solve the discrete form of the eigenproblem in Def. 5, i.e.

[TABLE]

where $\boldsymbol{A}_{j}^{o}$ is the discretisation matrix assembled in the overlap region and $\boldsymbol{X}_{j}$ is the discrete form of the partition of unity.

The matrix $\boldsymbol{A}_{j}$ has to be assembled with Dirichlet constraints on the domain boundary as prescribed by the given PDE problem. However, in contrast to the matrices needed for the one-level component of the two-level additive Schwarz method, no Dirichlet constraints are imposed on subdomain boundaries.

For assembling $\boldsymbol{A}_{j}^{o}$, the same boundary conditions can be applied. However, additionally, the matrix should only be assembled on the overlap region. Internally these elements are determined by adding a vector of ones across subdomains and checking for results greater than one.

The matrices $\boldsymbol{X}_{j}$ representing the partition of unity operator are diagonal and can be stored as vectors. Entries of $\boldsymbol{X}_{j}$ corresponding to Dirichlet domain boundaries or processor boundaries should be zero, and in sum they should add up to one across subdomains. Such a partition of unity is generated by adding vectors of ones with one communication between subdomains.

4.4 Solving the Eigenproblem

As the eigenproblems are defined per-subdomain, the eigensolver itself does not need to run in parallel. However, solving larger problems requires an efficient iterative solver. A suitable choice is ARPACK [8].As the eigenvalues of interest for the GenEO coarse space are those of smallest magnitude, the Shift and Invert Spectral Transformation Mode supported by ARPACK is used. Instead of the generalized eigenproblem $\boldsymbol{A}x=\boldsymbol{M}x\lambda$, ARPACK solves the transformed problem $(A-\sigma M)^{-1}Mx=x\nu$. The eigenvalues of the transformed problem are related to those of the original problem by $\nu=\frac{1}{\lambda-\sigma}$ and the eigenvectors are identical. In the transformed problem, the eigenvalues of the original problem whose absolute values are closest to $\sigma$ are now the eigenvalues of largest magnitude, and can therefore be efficiently solved by the Krylov method. Choosing $\sigma$ near zero, the method delivers the eigenvalues of smallest magnitude at good performance. Finally, in order to form the actual basis vectors, the eigenvectors are multiplied by $\boldsymbol{X}_{j}$ and then normalized in the $l^{2}$ norm, as ARPACK delivers vectors of strongly varying norms.

4.5 Scalable Coarse Setup

Assembling the coarse matrix $\boldsymbol{R}_{H}\boldsymbol{A}_{H}\boldsymbol{R}_{H}^{T}$ requires particular care, as it is a non-localized, not trivially scalable operation. Due to domain decomposition, the global matrix $\boldsymbol{A}$ is only available in distributed form as matrices $\boldsymbol{A}_{j}$. Exploiting local support of basis functions, the coarse matrix $\boldsymbol{A}_{H}$ breaks down into

[TABLE]

We note that $\varphi_{i}\boldsymbol{A}_{i}\varphi_{j}$ is zero for $\Omega_{i}\cap\Omega_{j}=\varnothing$, leading to a sparse structure in $A_{H}$. Therefore, all rows $i$ of $\boldsymbol{A}_{H}$ associated to basis functions $\varphi_{i}$ can be computed on the associated process locally while only requiring basis functions $\varphi_{j}$ from adjacent subdomains. In the implementationmultiple basis functions are communicated in a single step.

The resulting blocks are combined into a matrix $\boldsymbol{A}_{H}$ available on all processes, using direct MPI calls, while exploiting sparsity. Communication effort obviously increases with the dimension of $V_{H}$. This is a direct consequence of how two-level preconditioners are designed, and a good balance between coarse space size and preconditioner performance must be found.

The restriction and prolongation operators $\boldsymbol{R}_{H}$ and $\boldsymbol{R}_{H}^{T}$ are also only available locally. In case of the restriction $\boldsymbol{R}_{H}v_{h}$ of a distributed vector $v_{h}\in V_{h,0}(\Omega)$, it holds

[TABLE]

Each row $i$ can be computed by the process associated to $\varphi_{i}$, and the rows can be exchanged among all processes via MPI_Allgatherv. Again, the communication effort increases with the dimension of $V_{H}$.

Finally, the prolongation $\boldsymbol{R}_{H}^{T}v_{H}$ of a global vector $v_{H}\in V_{H}$ fulfills

[TABLE]

Here, each part of the sum associated with a processor can be computed locally and combined by nearest-neighbor communication, scaling ideally.

5 Numerical Experiments

In this section we demonstrate the solvers salient features, including its high parallel scalability up to $15,360$ cores, its robustness to heterogeneous material parameterss and its applicability to different elliptic PDEs. With exception of the final large-scale experiment all numerical examples in this section have been computed using the Balena HPC cluster of the University of Bath. Balena consists of $192$ nodes each with two $8$-core Intel Xeon E5-2650v2 Ivybridge processors, each running at $2.6$ GHz and giving a total of $3072$ available cores.

5.1 GenEO Basis on Highly Structured Problems

With clearly structured problems, it can be visually seen that the GenEO coarse space systematically picks up inclusions or channels in the parameter distribution. In Figure 3 the coarse approximation error is shown for a Darcy problem on a square domain. Dirichlet conditions are set to one at the top and zero at the bottom, Neumann conditions are set at the remaining boundary and a high-contrast parameter distribution with jumps and channels as shown on the left. We see that each inclusion has an effect on the approximation error. Adding additional eigenvectors from each subdomain to the coarse basis removes some of those error sources, the next eigenvectors pick up the skyscrapers and with only 4 eigenvectors per subdomain most channels are picked up. A total of $16$ coarse basis functions is enough to almost entirely solve the given problem.

5.2 Demonstration of Robustness

Robustness with respect to parameter contrast can be demonstrated solving the same Darcy problem as in section 5.1. We choose a subdomain decomposition into 8 by 8 squares, a two-cell overlap region diameter and a total of 800 $Q_{1}$ elements in each direction. Figure 4 (left) shows the resulting condition number for increasing contrast when setting up a GenEO basis with various numbers of eigenvectors per subdomain. Clearly, the asymptotic robustness guaranteed by the analysis is achieved in practice.

When running the same setup with a parameter distribution of 40 horizontal equally thick layers, it becomes clear from Figure 4 (right) that robustness is achieved exactly at four eigenvectors per subdomain. That stems from the fact that four coarse basis functions (together with the contribution of the one-level Schwarz method) are sufficient to represent the five layers contained in that subdomain. Similar relations can be observed with other strongly structured parameter distributions as well.

5.3 Comparison to other solvers

In this section we compare the performance of various preconditioned CG solvers. For this test we consider a flat composite plate of size $[0,100\mbox{mm}]\times[0,20\mbox{mm}]$. The laminate is made up of $12$ composite layers stacked in a sequence of different angles, refered to as a stacking sequence. The composite layers are seperated by very thin layers of resin. There is a large jump in material strength between the composite and resin layer and due to the rotated layers the anisotropy cannot be grid aligned. We discretise with quadratic, 20-node serendipity elements to avoid shear locking and use full Gaussian integration.

In Table 1 we compare the convergence of several iterative solvers. We record the condition number, the dimension of the coarse space $\mbox{dim}(V_{H})$ if applicable and the number of CG iterations required to achieve a residual reduction of $10^{-5}$. As expected the iteration counts increase steadily with the number of subdomains when no coarse space is used. In contrast, the iterations and the condition number estimates remain constant for the GenEO preconditioner as predicted by Thm. 3.1.

To demonstrate the efficiency of this approach we also compared with two different implementations of AMG. The implementation included in dune-istl has not originally been designed for composite application, and in its current form it does not seem to be robust, especially in parallel [9]. We applied this solver in its block version and measured strength of connection between two blocks for the aggregation procedure in the Frobenius norm. As a smoother we applied two SSOR iterations. In the test setup used here the dune-istl AMG converges very slowly or not at all, thus we do not include it in Table 1. We include results for BoomerAMG [15]. Here we retained the defaults for most parameters (HMIS coarsening without aggressive refinement levels and a hybrid Gauss-Seidel smoother). As recommended for elasticity problems we used blocked aggregation with block size corresponding to the spatial dimension. A strong threshold of $0.75$ was chosen after testing values in the range from $0.4$ to $0.9$. Due to a lower setup cost the boomerAMG solver is faster in actual CPU time than the GenEO solver in this small test case. However, for more complex geometries, boomerAMG does not perform very well and in our tests it does not scale beyond about $100$ cores in composite applications [7].

5.4 Industrially-motivated Example: Wingbox

In this section we describe an industrially motivated example in which we asses the strength of an airplane wingbox with a small localised wrinkle defect in one corner. Wrinkle defects often form during the manufacturing process and lead to strong local stress concentrations, which may cause premature failure [9, 10]. More details on this test setup can be found in [7].

We model a single bay of a wingbox of width $W=1$m, height $H=300$cm, length $L=1$m and internal radius of the corners $15$mm, as shown in Figure 5 (left). As in a typical aerospace application, the stacking sequence differs in the covers (top and bottom), corners and in the spar (sides). The changing stacking sequence is shown in Figure 5 (right). In total the laminate is made up of $39$ fibrous layers and $38$ resin layers. One of the corner radii contains a localised wrinkle with a parametrisation matching an observed defect in a CT-Scan of a real corner section, for further details see [10].

Two forms of loading are applied, an internal pressure of $0.109$MPa, arising from the fuel, is applied to the internal surface and a thermal pre-stress induced by the manufacturing process is imposed. The influence of ribs that constrain the wingbox in the $y$ direction are approximated by clamping all degrees of freedom at one end and tying all other degrees of freedom at the other end using a multipoint constraint.

To demonstrate the effectiveness of the GenEO solver we carry out a weak scaling and a strong scaling experiment. The experiments in this section were performed using the UK national HPC cluster Archer, which consists of $4,920$ Cray XC30 nodes with two $2.7$ GHz, $12$-core E5-2697 v2 CPUs each.

For the weak scaling experiment we refine the mesh, doubling the number of elements as we double the number of cores. Tab.  2 (left) contains the number of degrees of freedom, iteration numbers for the preconditioned CG, an estimate of the condition number of the discretisation matrix $\kappa(\boldsymbol{A})$, the dimension of the coarse space $\dim V_{H}$, as well as the total run time for each test. The weak scaling of the iterative CG solver with GenEO preconditioner is indeed almost optimal up to at least $15,360$ cores. We note that as shown in Thm. 3.1 the condition number remains bounded. As expected we observe a clear connection between the number of iterations and the condition number.

Tab.  2 (right) shows a strong-scaling experiment. The iterative CG solver with GenEO preconditioner scales almost optimally to at least $11,320$ cores. Memory constraints prevented tests with fewer than $2880$ cores. Tab. 2 shows that the number of iterations indeed remains almost constant. The last column gives the parallel efficiency, it remains high up to $11,320$ cores.

6 Acknowledgements

This work was supported by an EPSRC Maths for Manufacturing grant (EP/K031368/1). This research made use of the Balena High Performance Computing Service at the University of Bath. This work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk).