An Effective ${\mathcal{H}}^2$-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation P

Abstract

We present iterative solution strategies for solving efficiently Magneto-Quasi-Static (MQS) problems expressed in terms of an integral formulation based on the electric vector potential. Integral formulations give rise to discrete models characterized by linear systems with dense coefficient matrices. Iterative Krylov subspace methods combined with fast compression techniques for the matrix-vector product operation are the only viable approach for treating large scale problems, such as those considered in this study. We propose a fully algebraic preconditioning technique built upon the theory of H2-matrix representation that can be applied to different integral operators and to changes in the geometry, only by tuning a few parameters. Numerical experiments show that the proposed methodology performs much better than the existing one in terms of ability to reduce the number of iterations of a Krylov subspace method, especially for fast transient analysis.

1. Introduction

In this work, we consider the numerical solution of the Magneto-Quasi-Static (MQS) problem expressed in terms of an integral formulation based on the electric vector potential [1, 2, 3]. Integral formulations give rise to discrete models characterized by linear systems with dense coefficient matrices. Iterative Krylov subspace methods combined with fast compression techniques for the matrix-vector product operation are the only viable approach for treating large scale problems, such as those considered in this work. In this framework, a key role is played by the choice of the preconditioner, which has to be (i) robust enough to reduce significantly the number of iterations required to converge toward an accurate numerical solution, and (ii) efficient to reduce the overall computational and memory costs for the iterative solution. In this work, we focus on this critical computational aspect for the MQS problem. We propose a fully algebraic preconditioning technique built upon the theory of ${\mathcal{H}}^2$-matrix representation, which can be applied to different integral operators and to changes in the geometry only by tuning a few parameters.

Many algebraic preconditioning methods for solving dense linear systems arising from the discretization of integral operators are computed from a sparse approximation of the stiffness matrix, which retains the near-field mesh interactions giving rise to the most relevant contributions to the singular integrals of the model and is easy to be factorized or inverted [4, 5, 6]. However, our past experience in using such methods on our MQS problem, e.g., multilevel incomplete LU factorizations [7], sparse approximate inverses [8], and preconditioners computed from the sparse resistance matrix of the discrete system of equations (see Eq. 5 below) was not completely satisfactory. We found that the lack of global information due to the compact support of the sparse near-field matrix often results in slow convergence. For instance, in the CARIDDI code that implements the electric vector potential formulation of our problem [1, 9], a preconditioner based on the sparsity pattern arising from local interactions among neighboring basis functions is used to accelerate the iterative solution (see [10]). The aforementioned limitations are clearly shown in our numerical experiments, suggesting the need to develop global preconditioners that can incorporate contributions from the “far-field” interactions for this class of problems. In our model, far-field interactions can be included by taking into account the vector potential.

Prompted by these considerations, in this paper, we propose a class of preconditioners based on the so-called Hierarchical or $\mathcal{H}$-matrices [11], which provide a general mathematical framework for a highly compact and kernel independent accurate representation of integral equations with almost linear complexity [12]; see [13] for an example of application. $\mathcal{H}$ -matrices can be used to define global preconditioners built by means of the decomposition of small rank blocks at a moderate memory requirement. In previous works by other authors, highly efficient approximate LU-factorization and multilevel sparse approximate inverses preconditioners of boundary element matrices requiring an almost linear computational cost have been developed (see [14]) by means of the fast arithmetic of $\mathcal{H}$ -matrices. The novelty of this work is to use the latest development of $\mathcal{H}$ -matrices, that is, ${\mathcal{H}}^2$ -matrix theory [15, 16].

The paper is organized as follows. In Section 2, we present the MQS volume integral formulation used in our study and its numerical implementation. In Section 3, we establish the need for an effective preconditioner for our problem. In Section 4, we outline some fundamental properties of ${\mathcal{H}}^2$-matrix theory, and in Section 5, we describe the computational steps of a matrix factorization that has $\mathcal{O}\left(N\cdot log\left(N\right)\right)$ complexity, $N$ being the number of DoFs (Degrees of Freedom) and is applied for the design of our preconditioner. We assess the effectiveness of the proposed method for eddy current computations in fusion devices simulations by numerical experiments in Section 6, also against other algebraic preconditioning strategies. Finally, in Section 7, we draw some concluding remarks arising from the study and some perspective for future work.

2. The volume integral numerical formulation

Here we recall the main assumptions at the basis of a Magneto-Quasi-Static volume integral formulation and of its numerical implementation. We refer to a conducting region $V_C$, without any topological constraint, characterized by the ohmic constitutive equation $\mathbf{J}=\sigma \mathbf{E}$, where $\mathbf{J}$ is the current density, $\sigma$ is the material conductivity, and $\mathbf{E}$ is the electric field induced in the conductor by a set of time-varying external sources distributed with assigned current density ${\mathbf{J}}_S$.

The relevant equations in the magneto-quasi-stationary limit are

$$\mathbf{E}=-\frac{\partial \mathbf{A}}{\partial t}-\nabla \phi$$E1

$$\nabla \times \mathbf{A}=\mathbf{B}$$E2

where

$$\mathbf{A}\left(\mathbf{r},t\right)=\frac{{\mu }_0}{4\pi }{\int }_{V_c}\frac{\mathbf{J}\left({\mathbf{r}}^{\text{'}},t\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\text{'}}\right|}d{V}^{\prime }+{\mathbf{A}}_s\left(\mathbf{r},t\right)$$E3

and

$${\mathbf{A}}_s\left(\mathbf{r},t\right)=\frac{{\mu }_0}{4\pi }{\int }_{{\mathrm{ℝ}}^3-V_c}\frac{{\mathbf{J}}_s\left({\mathbf{r}}^{\text{'}},t\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\text{'}}\right|}d{V}^{\prime }.$$E4

By imposing the Ohm law in weak form, we have

$$\underset{V_c}{\int }\mathbf{W}\cdot \left\{{\sigma }^{-1}\mathbf{J}+\frac{\partial }{\partial t}\left[\frac{{\mu }_0}{4\pi }{\int }_{V_c}\frac{\mathbf{J}\left({\mathbf{r}}^{\text{'}},t\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\text{'}}\right|}d{V}^{\prime }\right]+\frac{\partial {\mathbf{A}}_s}{\partial t}\right\}\mathit{dV}=0$$E5

for any $\mathbf{W}\in S$ and with $\mathbf{J}\in S$, $S$ being the set of solenoidal vector functions with continuous normal component across $V_C$ and zero normal component on the boundary of $V_C$.

The numerical model implemented in the CARIDDI code [1, 9] is obtained by applying the Galerkin method to Eq. (5), having expanded the unknown currents $\mathbf{J}$ as $\mathbf{J}\left(\mathbf{r},t\right)={\sum }_k{I}_k\left(t\right)\nabla \times {\mathbf{N}}_k\left(\mathbf{r}\right)$, where ${\mathbf{N}}_k$ is the set of edge element basis functions whose uniqueness is assured by the tree-cotree gauge [17]. Upon discretization, we obtain the following linear dynamical system:

$$\left\{\begin{array}{l}\mathbf{L}\frac{d\mathbf{I}}{\mathit{dt}}+\mathbf{RI}=-\frac{d{\mathbf{V}}_0}{\mathit{dt}},t\ge 0\\ \mathbf{I}\left(0\right)={\mathbf{i}}_0\end{array}\right.,$$E6

where $\mathbf{I}\left(t\right)=\left[I_i\left(t\right)\right],\mathbf{R}=\left[R_{\mathit{ij}}\right],\mathbf{L}=\left[L_{\mathit{ij}}\right],{\mathbf{V}}_0\left(t\right)=\left[V_{0,i}\left(t\right)\right]$ and

$$R_{\mathit{ij}}=\underset{V_c}{\int }\nabla \times {\mathbf{N}}_i\left(\mathbf{r}\right)\cdot {\sigma }^{-1}\nabla \times {\mathbf{N}}_j\left(\mathbf{r}\right)\mathit{dV},$$E7

$$L_{\mathit{ij}}=\frac{{\mu }_0}{4\pi }\underset{V_c}{\int }\underset{V_c}{\int }\frac{\nabla \times {\mathbf{N}}_i\left(\mathbf{r}\right)\cdot \nabla \times {\mathbf{N}}_j\left({\mathbf{r}}^{\text{'}}\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\text{'}}\right|}\mathit{dVd}{V}^{\prime },$$E8

$$V_{0,k}\left(t\right)=\underset{V_c}{\int }\nabla \times {\mathbf{N}}_k\left(\mathbf{r}\right)\cdot {\mathbf{A}}_s\left(\mathbf{r},t\right)\mathit{dV}.$$E9

The quantity ${\mathbf{i}}_0$ in (6) denotes the prescribed initial condition. The vector potential ${\mathbf{A}}_s\left(\mathbf{r},t\right)$ is produced by a known source current density ${\mathbf{J}}_s\left(\mathbf{r},t\right)$_, and it can be computed by means of the Biot-Savart law. For the sake of simplicity, we assume that the source current density ${\mathbf{J}}_s$ is due to the current $i_s\left(t\right)$ flowing in a given coil.

The ODE in (6) is integrated in time by the following scheme [18]:

$$\left\{\begin{array}{l}{\mathbf{AI}}_n={\mathbf{LI}}_{n-1}+{\mathbf{V}}_{0,n}-{\mathbf{V}}_{0,n-1}\\ {\mathbf{I}}_0={\mathbf{i}}_0\end{array}\right.,$$E10

where ${\mathbf{I}}_n=\mathbf{I}\left(n\mathrm{\Delta }t\right)$ and ${\mathbf{V}}_n=\mathbf{V}\left(n\mathrm{\Delta }t\right)$, and matrix $\mathbf{A}$ is expressed as

$$\mathbf{A}=\mathbf{R}\mathrm{\Delta t}+\mathbf{L}.$$E11

The quantity $\mathrm{\Delta }t$ is the time integration step. Matrix $\mathbf{R}$ is sparse and local, whereas $\mathbf{L}$ is fully populated because it arises from the discretization of an integral operator. Note that $\mathbf{R}$, $\mathbf{L}$_, and $\mathbf{A}$ are positive definite. The solution of large linear systems with a dense coefficient matrix poses a relevant computational challenge as it requires fast matrix solvers with reduced algorithmic and memory complexity. In this work, we solve Eq. (10) iteratively, for each prescribed $n$. To speed up the iterative solution, we develop an accurate and computationally efficient preconditioner $\mathbf{P}$ built upon the ${\mathcal{H}}^2$ -matrix representation of matrix $\mathbf{A}$.

3. The need of an effective preconditioner

As it is well known, the main operations that are performed at each step of an iterative solver are: (i) the computation of the matrix-by-vector product $\mathbf{Ax}$, and (ii) the application of the preconditioner. The overall computational cost of the iterative solution is then proportional to the number of arithmetic operations required to carry out the multiplication $\mathbf{Ax}$ and to the number of iterations that are necessary to converge to a user-defined accuracy. The role of the preconditioner is to reduce the number of iterations and, ultimately, the overall solution cost [19]. In this work, we apply the preconditioner from the left and replace the linear system in (10) with the new equivalent system

$${\mathbf{P}}^{-1}{\mathbf{AI}}_n={\mathbf{P}}^{-1}\left({\mathbf{LI}}_{n-1}+{\mathbf{V}}_{0,n}-{\mathbf{V}}_{0,n-1}\right),$$E12

where the preconditioned stiffness matrix ${\mathbf{P}}^{-\mathbf{1}}\mathbf{A}$ is assumed to be well conditioned and to have most of its eigenvalues clustered around point one of the spectrum. Therefore, the “ideal” preconditioner has to be a good approximation to ${\mathbf{A}}^{-1}$. A critical issue in this context is the dependence of matrix $\mathbf{A}$ from $\mathrm{\Delta }t$. In fast transient analysis, $\mathrm{\Delta }t$ has to be small enough to model properly the underlying dynamics. On the other hand, in slow transients $\mathrm{\Delta }t$ has to be large enough to cover the full-time interval of interest at acceptable computational cost. The dependence of $\mathbf{A}$ on $\mathrm{\Delta }t$ yields the dependence of $\mathbf{P}$ on $\mathrm{\Delta }t$. According to (8), for small $\mathrm{\Delta }t$ the preconditioner should be tailored mainly on $\mathbf{L}$, whereas for large $\mathrm{\Delta }t$ the preconditioner should be tailored mostly on $\mathbf{R}$. Here and in the following, “small” and “large” $\mathrm{\Delta }t$ should be intended as compared with the slowest eigenmodes of the dynamical matrix ${\mathbf{L}}^{-1}\mathbf{R}$. In intermediate cases ($\mathrm{\Delta }t$ neither too small nor too large), the preconditioner should depend on both $\mathbf{R}$ and $\mathbf{L}$. Since the use of an implicit time stepping does not impose specific restrictions on $\mathrm{\Delta }t$, the appropriate value of $\mathrm{\Delta }t$ depends on the underlying dynamics induced by the source, i.e., on the rate of change of ${\mathbf{V}}_0\left(t\right)$ defined in (9) and appearing in (6).

A canonical waveform for $i_s\left(t\right)$, which can be used as the elementary building block to approximate more complex waveforms¹, is given by a linear ramp $i_s\left(t\right)=\mathit{\alpha t}$, for $t\ge 0$, being $\alpha$ the rate of change for the electrical current. Consequently, ${\mathbf{V}}_0\left(t\right)=\mathbf{W}t$ for $t\ge 0$ is also linear. The source term ${\mathbf{V}}_{0,n}-{\mathbf{V}}_{0,n-1}$ is equal to $\mathbf{W}\mathrm{\Delta }t$, which is constant for any $n$. Let us consider a transient starting from a vanishing initial condition ${\mathbf{i}}_0=\mathbf{0}$ in (6). As long as the current at the previous step is negligible, Eq. (10) reduces to ${\mathbf{AI}}_n=\mathbf{W}$. Therefore, the matrix to be inverted is $\mathbf{A}$_, and the preconditioner $\mathbf{P}$ has to be “tuned” on such matrix. On the other hand, when the stationary state is achieved, the solution ${\mathbf{I}}_n$ is constant, and Eq. (10) reduces to $\left(\mathrm{\Delta }t\mathbf{R}+\mathbf{L}\right){\mathbf{I}}_{\infty }={\mathbf{LI}}_{\infty }+\mathbf{W}$, that is, $\mathrm{\Delta }t{\mathbf{RI}}_{\infty }=\mathbf{W}$. In this case, the matrix to be inverted is $\mathbf{R}$_, and the preconditioner needs to be “tuned” on $\mathbf{R}$. This explains the reason why the preconditioner based on matrix $\mathbf{R}$ is not effective in dealing with fast dynamics. As an example, we measured its performances on a 3D transient eddy current problem where the driving current has a time-varying waveform given by $i_s\left(t\right)=\mathit{\alpha t}$, for $t\ge 0$, where $\alpha ={10}^7\mathit{Hz}$. This electrical current flows in a circular ring ($R=10\mathit{cm}$), centered on the axis of a flat conducting disk (radius $R_0=1m$, thickness $3\mathit{mm}$). Validation is carried out against the Femlab$\copyright$ code [Version 3.51], as shown in Figure 1 (right). Figure 1(left) shows the induced current distribution at $t=3.4\mathit{ms}$ together with the finite element mesh used in the computation and the ohmic losses in the conductor. Figure 2 shows the characteristics of the iterative inversion using $\mathbf{P}=\mathbf{R}$ as a preconditioner and the Generalized Minimal Residual Method (GMRES) [20] as the iterative solver. Since the number of discrete unknowns for this test problem is 22,657, the results of Figure 2 clearly show that that the resistive preconditioner $\left(\mathbf{P}=\mathbf{R}\right)$ is not effective at the early stage of the transient phase, when the inductive effects accounted by matrix $\mathbf{L}$ are not negligible.

Figure 1.

Left: discretization of the plate used for assessing the performances of the preconditioner, together with the induced current distribution due to a current on a coaxial ring in a significant instant of the transient. Right: power losses on the plate as a function of time computed by CARIDDI (blue) and by Femlab$\copyright$ code (red), under the hypothesis of axisymmetry.

Figure 2.

Number of iterations of GMRES using $R$ as a preconditioner during the transient.

4. ${\mathcal{H}}^2$-matrix representation

The ${\mathcal{H}}^2$-matrix representation is widely used for storing efficiently a dense matrix arising from the discretization of a boundary integral equation [15]. In ${\mathcal{H}}^2$-matrix theory, the degrees of freedom (DoFs) of the underlying finite element mesh are partitioned into small subsets (or cluster) of nodes (they are edges in our problem). We illustrate the partitioning process with reference to the geometry shown in Figure 3 that will be described in detail in the next section. At the first step, the set of DoFs is split in two clusters of nodes (see Figure 3 (left)); then, each cluster is split recursively (see Figure 3 (center) and (right)) until a minimum number of DoFs (called leaf size or LS) is reached by the recursive partitioning. The binary split is carried out using geometric information from the mesh. First, the coordinates of the barycenter of the DoFs are computed along the $x$, $y$_, and $z$ axes. Then, the axis (either $x$ or $y$ or $z$) corresponding to the largest geometric coordinate in magnitude is selected, and the DoFs are split accordingly by introducing a plane orthogonal to this axis and passing through the median of the distribution of the coordinates of the DoFs along the same axis (see Figure 4). A key point of the partitioning is that cluster of nodes are built to contain DoFs that are geometrically neighbors in the finite element mesh. From a geometric viewpoint, all the DoFs belonging to a given cluster of edges are contained within a parallelepiped, which is parallel to the coordinate planes (see Figure 5). The key parameters related to a cluster $N_k$ are: (i) the number of DoFs contained in the cluster, (ii) the parallelepiped $P_k$ containing the DoFs of $N_k$, (iii) the diameter $d_k=\left|P_k\right|$ and the center $c_k$ of $P_k$.

Figure 3.

The hierarchical partitioning of the DoFs. Left: The DoFs of the mesh are split into two sets (or clusters) of nodes. Center: The DoFs of the red set at the first step are partitioned in two subsets. Right: The DoFs of the red set at the second step are further partitioned into two subsets and so on.

Figure 4.

Example of DOFs splitting for a two-dimensional (2D) unstructured mesh.

Figure 5.

Top: Two clusters that are classified as far. The size of the individual regions (given by the diagonal of the parallelepiped that bounds the region) occupied by the DoFs is larger than the distance between the barycenter of the clusters. Bottom: Two clusters that are classified as near. The size of the individual regions occupied by the DoFs is not large enough if compared with the distance between the barycenter of the clusters.

This partitioning process can be efficiently represented in a binary tree data structure. The binary tree is of paramount importance to classify interactions between clusters of DoFs in terms of near (or inadmissible) interactions and far (or admissible) interactions. An interaction is defined far if the two clusters of DoFs are geometrically well separated in the mesh, and near otherwise, as shown in Figure 5. We will use a simple criterion based on geometric distance for the classification. Two cluster of nodes $N_k$ and $N_j$ are named far if the ratio

$${\xi }_{\mathit{kj}}\widehat{=}\frac{\left|{\mathbf{c}}_k-{\mathbf{c}}_j\right|}{max\left(d_k,d_j\right)},$$E13

is greater than a suitable threshold ${\xi }_{\mathit{THR}}>1.$

The mesh partitioning into clusters of nodes is represented by means of the so-called Block Cluster Tree (BCT) data structure (see Figure 6 and [16] for details). Thanks to its binary structure, the BCT can be easily mapped into a corresponding partitioning of the discretization matrix $\mathbf{L}$ and, hence, of $\mathbf{A}$ into small blocks, yielding naturally a decomposition $\mathbf{A}={\mathbf{A}}_{\mathit{\text{near}}}+{\mathbf{A}}_{\mathit{far}}$.

Figure 6.

Example of cluster tree partitioning ($n=8$ clusters). Left: The block cluster tree. Right: Matrix representation of the block cluster tree after four levels of matrix block partitioning. Admissible blocks are green, and inadmissible blocks are red.

Interactions of near clusters lead to high or full-rank blocks (the near part ${\mathbf{A}}_{\mathit{\text{near}}}$ of $\mathbf{A}$) that need to be explicitely computed and stored, whereas far clusters interactions give rise to low-rank blocks (the far part ${\mathbf{A}}_{\mathit{far}}$ of $\mathbf{A}$). In the ${\mathcal{H}}^2$ -matrix representation, the far interactions between, say, the $i$-th and $j$-th clusters, are approximated by a low-rank matrix factorization of the form

$${\mathbf{A}}_{\mathit{ij}}\approx {\mathbf{V}}_i\cdot {\mathbf{S}}_{\mathit{ij}}\cdot {{\mathbf{V}}_j}^{\mathrm{T}}.$$E14

Such factorization can be obtained by using a degenerate Lagrange polynomial approximation of the kernel [3]. This is possible because the kernel defined in Eq. (8) is smooth when it is evaluated on the far field. The columns of matrix ${\mathbf{V}}_i$ represent the expansion basis for the DoFs associated to the $i$-th cluster, and ${\mathbf{S}}_{\mathit{ij}}$ is the coupling matrix of dimensions $k_i\times k_j$ where $k_i$ and $k_j$ are the column dimensions of ${\mathbf{V}}_i$ and ${\mathbf{V}}_j$, respectively. It is worth noting that for far interactions, integers $k_i$ and $k_j$ are expected to be much smaller than the number of rows of matrices ${\mathbf{V}}_i$ and ${\mathbf{V}}_j$. This representation of the far field is highly efficient in terms of both memory and computational costs, since it involves only small-rank matrices (${\mathbf{V}}_i$ for any $i$, and ${\mathbf{S}}_{\mathit{ij}}$ for any $i$ and $j$). Indeed, linear memory complexity can be achieved by means of the ${\mathcal{H}}^2$ -matrix representation [12].

Hereafter, the following notation will be used: we denote by $N_p$ the number of Lagrange points per direction used for the degenerate kernel approximation, and by $e_{\mathit{ij}}$ the relative error of the low-rank representation of a matrix block ${\mathbf{A}}_{\mathit{ij}}$ associated to the interactions of the $i$-th and $j$-th far clusters:

$$e_{\mathit{ij}}=\frac{‖{\mathbf{A}}_{\mathit{ij}}-{\mathbf{V}}_i{\mathbf{S}}_{\mathit{ij}}{{\mathbf{V}}_j}^{\mathrm{T}}‖}{‖{\mathbf{A}}_{\mathit{ij}}‖},$$E15

where $\parallel \cdot \parallel$ is the conventional 2-norm for matrices.

5. The ${\mathcal{H}}^2$-matrix factorization

To develop an effective preconditioner $\mathbf{P}$ for solving Eq. (12), first we set $\mathbf{P}$ equal to ${\mathbf{A}}_{{\mathcal{H}}^2}$, the ${\mathcal{H}}^2$-matrix representation of $\mathbf{A}$. Then, by exploiting the hierarchical structure of ${\mathbf{A}}_{{\mathcal{H}}^2}$, we build an efficient LU factorization of $\mathbf{P}$. A core operation of the factorization of $\mathbf{P}$ is the computation of the Shur complement for the block matrix ${\mathbf{A}}_{{\mathcal{H}}^2}$. In general, this computation is very expensive. However, the ${\mathcal{H}}^2$-approximation introduces many off-diagonal blocks, which are vanishing upon a suitable orthonormal transformation of ${\mathbf{A}}_{{\mathcal{H}}^2}$. Hence, the overall factorization costs of $\mathbf{P}$ can be significantly reduced to $\mathcal{O}\left(N\cdot log\left(N\right)\right)$, $N$ being the number of DoFs [16, 21].

The factorization algorithm proceeds level by level, starting from the last level of the BCT and computing a partial LU decomposition at each cluster of the tree. With reference to the leaf level, the last level in the BCT, the original stiffness matrix $\mathbf{P}$ is expressed as:

$$\mathbf{P}={\mathbf{L}}_{\mathbf{o}}{\mathbf{C}}_1{\mathbf{U}}_{\mathbf{p}}$$E16

where ${\mathbf{C}}_1=\mathit{diag}\left({\mathbf{I}}_{d_1},{\mathbf{I}}_{d_2},\dots ,{\mathbf{I}}_{d_{\#\left(\mathit{\text{leaf\hspace{0.17em}clusters}}\right)}},\mathbf{C}\right)$ and ${\mathbf{L}}_{\mathbf{o}}$ and ${\mathbf{U}}_{\mathbf{p}}$ can be easily inverted because they are lower and upper triangular matrices, respectively. In (16), ${\mathbf{I}}_{d_i}$ is an identity matrix of proper size, and $\mathbf{C}$ is a full matrix of dimension much smaller than the original matrix $\mathbf{P}$. After the factorization (16) has been computed at the current level, the algorithm continues at a higher level of the BCT by applying recursively the same process to matrix $\mathbf{C}$.

Two key computational aspects have to be carefully considered and monitored during the factorization. The first one is the accuracy of the approximation ${\mathbf{A}}_{{\mathcal{H}}^2}$ to $\mathbf{A}$, which is used to define the preconditioner $\mathbf{P}$. Obviously, the closer ${\mathbf{A}}_{{\mathcal{H}}^2}$ to $\mathbf{A}$, the more effective $\mathbf{P}$ is as a preconditioner for the iterative solver. The second aspect refers to the error introduced in the approximated factorization of $\mathbf{P}$. In the rest of the section, we discuss in detail the factorization of $\mathbf{P}$.

Specifically, matrix $\mathbf{P}$ is factorized in terms of a partial LU decomposition [16]. For doing this operation efficiently, using some properties of ${\mathcal{H}}^2$-matrix theory is of paramount importance to reduce the overall cost. The complete multilevel factorization algorithm is presented in Algorithm 1.

Algorithm 1: The iterative factorization algorithm.

1. Set ${\mathbf{A}}_{\mathit{\text{curr}}}={\mathbf{A}}_{{\mathcal{H}}^2}$.

2. for$k$=last level to the root do.

3.   for any cluster $i$ at level $k$do.

4.     Set ${\mathbf{A}}_{\mathit{\text{curr}}}={\mathbf{Q}}_i^H{\mathbf{A}}_{\mathit{\text{curr}}}{\mathbf{Q}}_i$, where ${\mathbf{Q}}_i$ is defined in Eq. (14).

5.     Compute the partial LU factorization of ${\mathbf{A}}_{\mathit{\text{curr}}}$ as.

   ${\mathbf{A}}_{\mathit{\text{curr}}}={\mathbf{L}}_i\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathbf{A}}}_{\mathit{\text{curr}}}\end{array}\right]{\mathbf{U}}_i$, where ${\mathbf{L}}_i$ and ${\mathbf{U}}_i$ are lower and upper triangular matrices, respectively.

6.     Set ${\mathbf{A}}_{\mathit{\text{curr}}}=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathbf{A}}}_{\mathit{\text{curr}}}\end{array}\right]$

7.   end for.

8. end for.

At step 5, matrix ${\mathbf{Q}}_i$ encompasses the basis related to the $i$-th cluster as follows:

$${\mathbf{Q}}_i=\mathit{diag}\left\{{\mathbf{I}}_{d_1},{\mathbf{I}}_{d_2},\dots ,\left[{\mathbf{V}}_i^{\perp },{\mathbf{V}}_i\right],\dots ,{\mathbf{I}}_{\#\left(\mathit{\text{leaf\hspace{0.17em}clusters}}\right)}\right\}.$$E17

The orthonormal transformation at step 5 is computationally cheap to carry out because of the simplified structure of matrix ${\mathbf{Q}}_i$, which involves only the $i$-th diagonal block of matrix ${\mathbf{A}}_{\mathit{\text{curr}}}$, of dimension $d_i\times d_i$. In the definition of ${\mathbf{Q}}_i$ in Eq. (17), ${\mathbf{V}}_i^{\perp }$ is the orthogonal complement of the cluster basis ${\mathbf{V}}_i$. This update of ${\mathbf{A}}_{\mathit{\text{curr}}}$, which essentially corresponds to a change of basis, is very effective on the far blocks to the $i$-th diagonal block of ${\mathbf{A}}_{\mathit{\text{curr}}}$. Thanks to (14), it is easy to verify that the first $d_i-k_{\mathrm{\ell }}$ rows of these blocks are zeroed out, where $k_{\mathrm{\ell }}$ is the rank used for the approximate factorization at the $\mathrm{\ell }$-th level of the BCT. Since most of the blocks in the original matrix $\mathbf{A}$ are far, they will be efficiently compressed by the procedure described above (see [16]). At step 10, a complete LU factorization of a dense matrix is required. However, this operation is rather cheap because matrix ${\tilde{\mathbf{A}}}_{\mathit{\text{curr}}}$, resulting from the previous steps, is small.

We highlight an important computational aspect of Algorithm 1. Upon processing the $i$-th cluster at step 5, the partial LU factorization of ${\mathbf{A}}_{\mathit{\text{curr}}}$ writes as

$${\mathbf{A}}_{\mathit{\text{curr}}}=\left[\begin{array}{cc}{\mathbf{A}}_{i^{\prime }{i}^{\prime }}& {\mathbf{A}}_{i^{\prime }p}\\ {\mathbf{A}}_{p{i}^{\prime }}& {\mathbf{A}}_{\mathit{pp}}\end{array}\right]=\underset{{\mathbf{L}}_i}{\underbrace{\left[\begin{array}{cc}{\mathbf{L}}_{i^{\prime }{i}^{\prime }}& \mathbf{0}\\ {\mathbf{A}}_{p{i}^{\prime }}{\mathbf{U}}_{i^{\prime }{i}^{\prime }}^{-1}& \mathbf{I}\end{array}\right]}}\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}& {\tilde{\mathbf{A}}}_{\mathit{\text{curr}}}\end{array}\right]\underset{{\mathbf{U}}_i}{\underbrace{\left[\begin{array}{cc}{\mathbf{U}}_{i^{\prime }{i}^{\prime }}& {\mathbf{L}}_{i^{\prime }{i}^{\prime }}^{-1}{\mathbf{A}}_{i^{\prime }p}\\ {\mathbf{A}}_{p{i}^{\prime }}& \mathbf{I}\end{array}\right]}}$$E18

where ${\mathbf{A}}_{i^{\prime }{i}^{\prime }}={\mathbf{L}}_{i^{\prime }{i}^{\prime }}{\mathbf{U}}_{i^{\prime }{i}^{\prime }}$, $i^{\prime }$ is equal to the number of nodes in the $i$-th clusters minus the approximation rank at the current level, and

$${\tilde{\mathbf{A}}}_{\mathit{\text{curr}}}={\mathbf{A}}_{\mathit{pp}}-{\mathbf{A}}_{p{i}^{\prime }}{\mathbf{U}}_{i^{\prime }{i}^{\prime }}^{-1}{\mathbf{L}}_{i^{\prime }{i}^{\prime }}^{-1}{\mathbf{A}}_{i^{\prime }p}.$$E19

Clearly, matrix ${\tilde{\mathbf{A}}}_{\mathit{\text{curr}}}$ still has a block structure, which is inherited from the cluster-to-cluster interactions. Contributions to the update term ${\mathbf{A}}_{p{i}^{\prime }}{\mathbf{U}}_{i^{\prime }{i}^{\prime }}^{-1}{\mathbf{L}}_{i^{\prime }{i}^{\prime }}^{-1}{\mathbf{A}}_{i^{\prime }p}$ due to near clusters mesh interactions are computed explicitely. On the other hand, contributions due to far clusters interactions are approximated using a low-rank factorization of the form given by Eq. (11), where cluster bases are simply updated by appending new columns pointing in the direction of the space orthogonal to the range of ${\mathbf{V}}_i$, thus avoiding to re-compute the whole basis from scratch. This approximate computation of the Schur complement is essential to achieve $\mathcal{O}\left(N\cdot log\left(N\right)\right)$ complexity. For further details, see [16].

Summing up, in this method two important sources of errors are introduced during the factorization:

• The far-field approximation. Rather than factorizing the original matrix $\mathbf{A}$, we factorize its ${\mathcal{H}}^2$-matrix representation ${\mathbf{A}}_{{\mathcal{H}}^2}$. The accuracy of this approximation depends on the number of Lagrange points ($N_p$) used to compute the ${\mathcal{H}}^2$-matrix approximation of $\mathbf{A}$.

• The error introduced in the Schur complement update (6) for each processed cluster. This error is controlled by a suitable threshold parameters and is the dominant one in the presented algorithm.

6. Numerical results and discussion

The numerical example described in this section refers to eddy current computations in fusion devices [22, 23]. The mesh of our case study is depicted in Figure 7. The plasma made of hydrogen isotopes inside the vacuum chamber of the device can be broadly described as a deformable current-carrying conductor in the frame of the so-called Magneto-Hydro-Dynamics (MHD) equations. Consequently, the electromagnetic interaction of the plasma with surrounding conducting structures is fundamental in the description of the behavior of such devices. In nominal situations, when one wants to intentionally control plasma current, position, and shape [24], this interaction occurs on a time scale of the order of tens or hundreds of milliseconds in the largest existing devices. Conversely, also off-normal events can occur, such as the so-called disruptions [25]; in such cases, the plasma energy content is released on a much shorter time (milliseconds or submilliseconds time range).

Figure 7.

The analyzed case study. Left: mesh of ITER magnet system. Right: A ${45}^{\circ }$ sector.

It should be noted that the geometry of fusion devices, which is ideally axisymmetric, in practice can show significant three-dimensional features (ports, holes, cuts, slits, etc.), which must be taken into account when computing the magnetic fields, the current densities, and the related force densities, which are required to get a satisfactory electromagnetic modeling. Summing up, fusion devices are particularly challenging from the computational electromagnetic viewpoint: they are multiphysics devices, with complicated three-dimensional features giving rise to huge computational models, with different timescales of interests involved. The methods developed in the present paper are designed to enable us to tackle such intrinsic difficulties to some extent. The number of nodes and elements of the finite element mesh is $249849$ and $132562$, respectively, resulting in $145243$ DoFs. The first three partitioning levels for such mesh are schematically depicted in Figure 3. In this study, the time step $\mathrm{\Delta }t$ is set equal to $1\mathit{ms}$ and $0.1\mathit{ms}$, to account for the different timescales mentioned above.

As aforementioned, the main operations required by an iterative solver are: (i) the computation of the matrix-by-vector product $\mathbf{Ax}$ and (ii) the application of the preconditioner at every iteration. We choose GMRES as the iterative solver. The preconditioner is the compressed approximate matrix ${\mathbf{A}}_{{\mathcal{H}}^2}$ presented in Section 5. Clearly, the use of a preconditioner introduces additional costs in terms of memory and computational cost, during the preprocessing and at each iteration. Therefore, it pays off only if there is a proper reduction of the number of required iterations. Here, to assess the performances of the preconditioner, we use matrix $\mathbf{A}$ without any approximation in computing the matrix-by-vector product during GMRES iterations. We stop the iterative process either when the initial preconditioned residual norm is reduced by five orders of magnitude or after 90 iterations. In our case, the computational time required to apply the preconditioner is proportional to its memory footprint. This is because of the almost linear complexity of the overall procedure (see [16] for details): both computational time and memory requirements increase linearly with the number of DoFs.

The performance of the proposed approach depends mainly on two key parameters: the number of Lagrange points $N_p$ and the leaf size $\mathit{LS}$. We remind that $\mathit{LS}$ determines the size of matrix ${\mathbf{A}}_{\mathit{\text{near}}}$ arising from near cluster interactions, whereas $N_p$ affects the number of columns of ${\mathbf{V}}_i$, the size of ${\mathbf{S}}_{\mathit{ij}}$ and, consequently, the accuracy of ${\mathbf{A}}_{\mathit{far}}$. It is worth noting that ${\mathbf{A}}_{\mathit{far}}$ in the ${\mathcal{H}}^2$ -matrix representation does not contribute significantly to the overall required memory: in all numerical experiments of this work, the memory required by ${\mathbf{A}}_{\mathit{far}}$ is less than 10% of the storage required by ${\mathbf{A}}_{\mathit{\text{near}}}$.

From (11), we have that the preconditioner memory footprint depends only on the storage of matrices ${\mathbf{L}}_0$, ${\mathbf{U}}_p$ and the dense but smaller block $\mathbf{C}$ of matrix ${\mathbf{C}}_1$ appearing in (16). Figure 8 shows the memory required by the preconditioner ($M_{\mathit{LU}}$) as a function of $\mathit{LS}$ and $N_p$. The result is normalized with respect to the memory required to store the full matrix $\mathbf{A}$, which is about $169\mathit{GB}$.

Figure 8.

Memory occupation as $\mathit{LS}$ and $N_p$ varies.

The memory footprint required by matrix ${\mathbf{C}}_1$ decreases with $\mathit{LS}$. Indeed, an increase of $\mathit{LS}$ enlarges the size of the block-to-block interactions, resulting in a smaller reduced matrix ${\mathbf{C}}_1$ because of a more effective compression of the block-to-block interactions. On the other hand, the sizes of ${\mathbf{L}}_0$ and ${\mathbf{U}}_p$ depend on: (i) the size of the near part, strictly depending on $\mathit{LS}$, and (ii) the ability of the ${\mathcal{H}}^2$ -matrix representation to compress the admissible block-to-block far-field interactions.

If we reduce $\mathit{LS}$, the memory required by the near interactions matrix ${\mathbf{A}}_{\mathit{\text{near}}}$ decreases, but the memory required by the far-field block-to-block interactions matrix increases. Indeed, we observe that at $\mathit{LS}=500$ the required memory depends on $N_p$, whereas for $\mathit{LS}=2000$ the required memory is almost independent on $N_p$. Eventually, we have an intermediate dependence for $\mathit{LS}=1000$. The drawback related to the memory demands can be addressed in a parallel environment.

To assess the effectiveness of the proposed preconditioner, we make a comparison with the “reference” preconditioner based solely on matrix $\mathbf{R}$. This choice resulted to be much more effective than many other state-of-the-art algebraic preconditioners, such as robust multilevel, incomplete LU factorizations, and sparse approximate inverse methods to accelerate the GMRES convergence. Also, we evaluated a preconditioner based on the sum of matrix $\mathbf{R}$ plus a contribution from matrix $\mathbf{L}$ on the same pattern of $\mathbf{R}$. This preconditioner failed to converge at increasing mesh size. Summing up, finding a proper preconditioner for this class of large-scale problems is challenging.

Figure 9 (top $\mathrm{\Delta }t=1\mathit{ms}$) and (bottom $\mathrm{\Delta }t=0.1\mathit{ms}$) shows the preconditioned relative residual norm $E_R$ of the GMRES solver versus the number of iterations for different $\mathit{LS}$ and $N_p$.

Figure 9.

Performance of ${\mathcal{H}}^2$ -LU preconditioner, for different values of the leaf size ($\mathit{LS}$) and number of Lagrange points ($N_p$). Top: $\mathrm{\Delta }t=1$ ms. bottom: $\mathrm{\Delta }t=0.1$ ms. performances from the reference preconditioner based on $R$ are also shown.

The numerical results clearly show that:

• The ${\mathcal{H}}^2$ -LU generally performs much better than the preconditioner based on $\mathbf{R}$ preconditioner.

• The efficiency of the preconditioner strongly depends on the choice of the parameters $N_p$ and $\mathit{LS}$ that control the accuracy of ${\mathbf{A}}_{\mathit{far}}$.

We highlight that the accuracy of ${\mathbf{A}}_{\mathit{far}}$ is rather low for $N_p=1$ and that the preconditioner based on $\mathbf{R}$ performs better in this situation. In the other cases, the proposed ${\mathcal{H}}^2$ -LU preconditioner outperforms the reference one based on $\mathbf{R}$, proving that it incorporates a good mechanism to balance contributions from matrices $\mathbf{L}$ and $\mathbf{R}\mathrm{\Delta }t$ in $\mathbf{A}$. This explains why at the smaller $\mathrm{\Delta }t$, when $\mathbf{L}$ gives the dominant contribution to $\mathbf{A}$, the number of iterations for the $\mathbf{R}$ based preconditioner increases. On the other hand, when using the ${\mathcal{H}}^2$ -LU preconditioner, which accounts for both $\mathbf{R}\mathrm{\Delta }t$ and the dominant part of $\mathbf{L}$, the number of iterations remains almost unchanged with $\mathrm{\Delta }t$.

The accuracy of ${\mathbf{A}}_{\mathit{far}}$ is a critical point for the effectiveness of the proposed preconditioner. To asses this, we define two quantities, denoted as $E_A$ and $E_{\mathit{SD}}$, that are the mean and the standard deviation for the relative error $e_{\mathit{ij}}$ defined in (15) arising from all the admissible block-to-block pairs $\left(i,j\right)$. As expected, Figure 10 shows that the error decreases as $N_p$ increases. Note that for $N_p=1$, a relative error larger than $0.1$ can be obtained. Eventually, the accuracy of the preconditioner is shown in Figure 11 for different values of $\mathit{LS}$, $N_p$_, and $\mathrm{\Delta }t$. To quantify the error, we introduce the metric $E_S$ defined as

Figure 10.

Mean and standard deviation of the relative error in $A_{\mathit{far}}$.

Figure 11.

Solution error $E_s$, as $\mathit{LS}$, $N_p$_, and $\mathrm{\Delta }t$ vary.

$$E_S=\frac{‖I_{\mathit{app}}-I‖}{‖I‖},$$E20

where $I$ is computed using a direct method, whereas $I_{\mathit{app}}$ is computed with the proposed preconditioned iterative solver. In (19), $‖\cdot ‖$ is the 2-norm, and both $I$ and $I_{\mathit{app}}$ are evaluated at a prescribed time.

Figure 11 confirms that the accuracy of the solution increases with $N_p$. Only for $\mathit{LS}=500$ and $\mathrm{\Delta }t=1\mathit{ms}$ the solution with $N_p=3$ is a little bit less accurate than the one with $N_p=2$. It is worth noting that as the accuracy of ${\mathbf{A}}_{\mathit{far}}$ increases, the iterative scheme converges in very few iterations. Summing up, these numerical results confirm that an efficient iterative solver for the CARIDDI integral model of the eddy current problem needs a preconditioner incorporating some far-field approximations from the interaction matrix $\mathbf{L}$.

7. Conclusions and future work

In this paper, a new preconditioner based on the ${\mathcal{H}}^2$ -LU decomposition has been introduced for solving integral equations models for eddy current problems. The application case refers to a large-scale model arising in nuclear fusion devices. The numerical results show that the preconditioner exhibits better performances compared with the purely resistive preconditioner. Specifically, the required number of iterations has been strongly reduced at different time steps. Moreover, the increase of costs in terms of both memory and CPU time during the preprocessing and at each iteration is largely paid off by the reduced number of iterations. In our view, two possible areas of future developments are: (i) deploying a parallel implementation of the ${\mathcal{H}}^2$ -LU decomposition in order to speed up the code execution and reduce the memory demands, and (ii) studying the performances of the method when it is used as a direct solver.

Acknowledgments

This work was partially supported by the program “Dipartimenti di Eccellenza 2018–2022,” MIUR, Italian Ministry of University and Research, by the project “CAR4PES,” 4th Marconi Fusion cycle, MUR PRIN Grant# 20177BZMAH and by the Research Südtirol/Alto Adige 2019 grant FH2ASTER, Provincia autonoma di Bolzano/Alto Adige – Ripartizione Innovazione, Ricerca, Università e Musei, contract nr. 19/34.