Introductory Chapter: Large Eddy Simulation for Turbulence Modeling

*Aamir Shahzad, Muhammad Kashif and Fazeelat Hanif*

## **1. Introduction**

Laminar and non-laminar are two major types of flows and are discussed mainly in fluid mechanics. Streamline and turbulent flows are examples of laminar and non-laminar flows; however, a laminar flow can be transformed into non-laminar by applying some kind of perturbations (such as temperature, pressure, force field, etc.) and/or by employing some gradient. This type of conversion process from laminar to non-laminar flow produces new patterns in between the two states, and these patterns are unstable, and it generates the flow instabilities in the fluid. Only large eddies (large-scale motions) are directly computed in large eddy simulation (LES); therefore, to create the three-dimensional (3D) unsteady governing equations for large-scale motions, a low-pass spatial filter is used for the instantaneous conservation equations. LES has wide range of applications for compressible flows, turbulent combustion, aeroacoustics, turbulent/transitional flows, and atmospheric sciences. In comparison to LES for cases involving incompressible flows, much less work has been performed on LES for compressible flows, and there are numerous difficulties/problems in this field. Extra work/requirements are needed for supersonic flows with shock waves in order to accurately and steadily capture the shock while also providing the spatial accuracy necessary to simulate a number of fine-scale turbulence structures. Low-order techniques are typically used to address shock waves, frequently using upwind schemes that are not particularly suitable for LES. Favre filtering is typically used in compressible flows to prevent the entry of sub-grid scale (SGS) terms into the continuity equation; therefore, knowledge and expertise obtained in incompressible flows may not be applicable. SGS modeling for compressible flows is significantly more difficult as a result of the additional equations that must be solved, such as the energy equation for the compressible case, and the necessity to represent additional SGS terms, such as the SGS heat flux [1].

With applications in a variety of combustion issues, LES of turbulent combustion first emerged in the 1990s and has grown significantly in the last 10 years. The majority of combustion chemistry takes place in SGS; therefore, models must be created because chemical reactions typically take place on sizes much smaller than those of LES meshes. However, even with very straightforward SGS combustion models, LES has showed considerable potential in this field and clearly outperformed the Reynolds-averaged Navier-Stokes (RANS) approach. However, due to the complexity of turbulent combustion, which includes chemical reactions, turbulence/chemistry interactions, liquid fuel atomization, liquid fuel injection, droplet breakup and evaporation, small-scale

molecular fuel air mixing, large-scale turbulent fuel air mixing, and chemical reactions in aircraft engines, there are significant challenges in this area. Many of these processes take place at various ranges of length and time [1]. LES also has wide range of applications in gas turbine. Reduced life cycle costs for the operator and decreased environmental effect during engine production and operation are crucial factors for gas turbine manufacturers. Many gas turbine components use a stable RANS approach; however, when it comes to the rotor-stator interaction in turbomachinery, we frequently use an unsteady RANS (URANS) approach where the governing flow equations are phase averaged assuming a constant rate of rotation [2]. LES-related techniques are being used more frequently as a result of the inability of RANS models to produce accurate off-design aerodynamics forecasts, noise source information, and predictive capability for the control of drag, noise, and mixing processes in general. The huge advances in accessible computing power over the past few years have provided additional encouragement for the usage of LES [3].

#### **1.1 LES for aeroacoustics**

A major amount of the noise that is produced by air and land transportation, such as fan noise, jet noise, high-speed train noise, and airframe noise, is a growing environmental concern. Turbulence is a prominent source of aerodynamic noise. Because LES directly computes large scale fluctuations, which are known to add the most to the noise generated in many issues, LES is a very helpful technique in aeroacoustics. Applications of LES for foretelling aerodynamic noise likely began in the 1990s and have since grown to be a very active area of research. LES shows considerable promise for aeroacoustics computations, from improving source modeling for acoustic analogies to practical prediction and designing of engineering systems in the near future. It also advances fundamental understanding of noise creation. LES algorithms should be able to reliably mimic the flow physics that captures the transfer of energy from turbulent to acoustic modes if properly built and validated. The proper SGS modeling, numerical issues such as high-order accuracy and careful usage of the boundary conditions, and practical engineering configurations where flow Reynolds numbers are typically very high make it impractical to use LES for both noise source capturing and its propagation and are all still significant challenges. Additionally, standard validation study against approved experimental databases can be carried out for relatively basic LES applications. Since intricate statistics such as two-point space-time correlations are essential to flow-generated sound, more attention should be given while validating LES applications in aeroacoustics, according to the theory. As a result, the validation may begin with the most basic facts before moving on to more intricate and acoustically important statistics [1].

## **1.2 LES for turbulent/transitional flows**

Turbulence has an important property, which is scale invariance. Certain characteristics of the flow are said to be scale invariant if they hold true across various motion scales. Such symmetry can be explained as an especially straightforward relation between small and large scales, making it a crucial component in turbulence models [4]. LES, a different strategy, was first suggested by Smagorinsky in 1963. The traditional RANS approach, which requires solving additional modeled transport equations to determine the so-called Reynolds stresses as a result of the averaging procedure, is not used by LES. When compared with direct numerical simulation (DNS), the computational cost of LES is significantly lower since only small-scale SGS

#### *Introductory Chapter: Large Eddy Simulation for Turbulence Modeling DOI: http://dx.doi.org/10.5772/intechopen.108294*

motions are represented and the large-scale motions of turbulent flow are computed directly. LES is more perfect than the RANS approach as the large eddies contain most of the turbulent energy and are accountable for most of the turbulent mixing and momentum transfer, and these eddies are captured directly by LES in full detail where they are modeled in the RANS approach. Moreover, the small scales tend to be more homogeneous and isotropic than the large ones, and thus, the SGS motions modeling should be easier than all scales modeling within a single model such as in the RANS approach [1]. There are two main types of SGS models for LES of turbulent flows. The models that yield expressions for SGS terminology such as a heat flux or stress tensor and typically include eddy viscosity notions fall under one group. The SGS stresses are secondary values, which are directly computed from the definitions in the other category, which describes the unresolved primitive variables such as velocity or temperature [5]. For most meteorological applications, large eddies, which hold the majority of the turbulent kinetic energy (TKE) also referred as energy-containing eddies, are the most significant scales for atmospheric planetary boundary layer (PBL) turbulence, as an example, and are in charge of most turbulent transport concerned. It is a simulation that specifically determines (or eliminates) large eddies when LES roughly reflects the effects of smaller ones. As LES's grid resolution increases, less are present, a greater spectrum of turbulent eddies is resolved, and LES-generated flows are parameterized, and they become more representative throughout the flow field. Thus, now LES is the most promising/feasible numerical mean for realistic turbulent/ transitional flows simulation [6]. LES inlet conditions were also generated on the basis of implementation of digital filter generator (DFG). The test case was the LES of a channel flow with a continuously repeated constriction. The DFG was used to construct three-time series that would be used as LES inlet conditions in the future. Along with entering the first and second moments of the velocity field over the inlet plane from the periodic boundary condition (PBC) simulation in the first, the DFG's input turbulence scales were set to be spatially homogeneous with values determined using a channel inlet height-weighted area average. The turbulence scales were allowed to change in the second and third time series. Their variation is again inferred from the PBC simulation and is related to wall normal direction. Then, LES's inlet boundary conditions were created using these distinct time series. These changing turbulent scales have increased the simulation accuracy, which has significant uses [7].

## **1.3 LES for atmospheric science**

Prior LES research mostly concentrated on a turbulent PBL with flat terrain and no clouds. Due to the enormous thermal plumes present and the lack of other flow conditions, this flow regime is best suited to LES without involving complex physical phenomena (such as latent heating and radiation). However, in recent years, LES has been broadened to study more challenging and complicated PBL regimes, which are pertinent to forecasts for severe weather and climate, or more recent applications for wind energy [6].

## **2. History, current state, and future challenges of LES**

## **2.1 History**

Smagorinsky first proposed LES in 1963 for forecasting air flow, and early uses likewise fell under this category. Deardoff and Schumann introduced LES to engineering-related flow for the first time in 1970 and 1975, respectively. From the 1960s to roughly the middle of the 1980s, LES took a while to develop, and the applications primarily consisted of straightforward, building-block flows, such as homogeneous turbulence, plane channel flows, mixing layers, and so on. However, as computing power increased, a very fast development and sharp increase in LES applications began around the middle of the 1980s, particularly after the 1990s with substantial growth of the LES community and a wide range of LES applications shifting from simple flows to complex flows, including multi-phase flow, heat transfer, and fluid dynamics, aeroacoustics, transmission, combustion, etc. In addition to the rise in computing power, it is now evident that RANS approaches naturally have limitations and can't deal with certain categories of sophisticated turbulent flow issues, and LES has developed quickly and has a wide range of applications [1].

## **2.2 Current state**

As was indicated in the preceding section, the LES was initially applied successfully to examine the specifics of flow problems with low Reynolds numbers and generally simple geometry, such as homogeneous turbulence, mixing layers, and flat channel flows. Despite the fact that LES is used in such academic or essential setups still exists today, primarily for model validation and a basic comprehension of flow physics, etc. When the RANS technique has failed, emphasis has changed to more intricate arrangements with flow characteristics. Particularly, corporate interest in applying LES to complicated engineering flows has been stimulated by decades of advancement in LES and the advent of cheap workstation clusters and massively parallel computers. Nevertheless, the RANS technique has not been substituted by LES, and it seems unlikely that it will be for some time due to two key factors, computational analytical tool for real-world engineering challenges firstly, notwithstanding the present computing ability for practical purposes, performing LES routinely still costs far too much in terms of computing; secondly, LES is not yet mature enough for users without sufficient results can be achieved with the level of solution accuracy that can be anticipated using experience and knowledge. Currently, LES has wide range of applications in turbulent flows, gas turbine, jet noise, aeroacoustics, and atmospheric science [1].

## **2.3 Future challenges**

In the future, LES will be properly used for a wider series of flow problems and for more difficult problems with more multi-disciplinary uses. Nevertheless, numerous major issues/challenges related with LES and its applications such as wall layer modeling, SGS modeling, LES of turbulent combustion, generation methods for inflow boundary conditions, etc., are still existing. Before LES can become a trustworthy, strong engineering analysis tool that can be utilized as a substitute for RANS, there are still important problems that need to be overcome. It is extremely improbable that LES will entirely overtake RANS and become a design tool for the foreseeable future, without considerable years of LES experience [1].

Actually, the PBL is more complex as compared with LES for simulation of systems; however, this complexity is generated due to heterogeneous nature of the considering surface. For instance, the earth land surface is described by spatially varying elements, urban expansion, and undulating grounds, which can merge circulations and therefore alter the turbulence dynamics. These complex surfaces may *Introductory Chapter: Large Eddy Simulation for Turbulence Modeling DOI: http://dx.doi.org/10.5772/intechopen.108294*

significantly affect turbulence transport in many climate applications, for example, vegetation development, pollution, cloud formation, and storm formation. In summary, the LES is more employed than realistic PBL for multi-dimensional environmental flows.

## **Abbreviations**


## **Author details**

Aamir Shahzad\*, Muhammad Kashif and Fazeelat Hanif Modeling and Simulation Laboratory, Department of Physics, Government College University Faisalabad (GCUF), Faisalabad, Pakistan

\*Address all correspondence to: aamirshahzad\_8@hotmail.com

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## **References**

[1] Zhiyin Y. Large-eddy simulation: Past, present and the future. Chinese Journal of Aeronautics. 2015;**28**(1):11-24

[2] Menzies. Large eddy simulation applications in gas turbines. Philosophical Transactions of the Royal Society A. 2009;**367**:2827-2838

[3] Tucker P. The LES model's role in jet noise. Progress in Aerospace Sciences. 2008;**44**(6):427-436

[4] Meneveau C, Katz J. Scale-invariance and turbulence models for large-eddy simulation. Annual Review of Fluid Mechanics. 2000;**32**(1):1-32

[5] Domaradzki JA, Adams NA. Direct modelling of subgrid scales of turbulence in large eddy simulations. Journal of Turbulence. 2002;**3**(1):024

[6] Moeng CH, Sullivan PP. Large-eddy simulation. Encyclopedia of Atmospheric Sciences. 2015;**2**:232-240

[7] Veloudis I, Yang Z, McGuirk JJ, Page GJ, Spencer A. Novel implementation and assessment of a digital filter based approach for the generation of LES inlet conditions. Flow, Turbulence and Combustion. 2007;**79**(1):1-24

## **Chapter 2**

## An Effective H<sup>2</sup> -LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation Problems

*Salvatore Ventre, Bruno Carpentieri, Gaspare Giovinco, Antonello Tamburrino, Fabio Villone and Guglielmo Rubinacci*

## **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 H<sup>2</sup> -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.

**Keywords:** H<sup>2</sup> -matrix, LU decomposition, MQS volume integral formulation, time domain, preconditioning, eddy currents, controlled thermonuclear fusion

## **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–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 H<sup>2</sup> -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–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 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. 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 H matrices. The novelty of this work is to use the latest development of H -matrices, that is, H<sup>2</sup> -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 H<sup>2</sup> -matrix theory, and in Section 5, we describe the computational steps of a matrix factorization that has Oð Þ *N* � log ð Þ *N* 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

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

region *VC*, without any topological constraint, characterized by the ohmic constitutive equation **J** ¼ *σ***E**, where **J** is the current density, *σ* is the material conductivity, and **E** is the electric field induced in the conductor by a set of time-varying external sources distributed with assigned current density **J***S*.

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

$$\mathbf{E} = -\frac{\partial \mathbf{A}}{\partial t} - \nabla \rho \tag{1}$$

$$\nabla \times \mathbf{A} = \mathbf{B} \tag{2}$$

where

$$\mathbf{A}(\mathbf{r},t) = \frac{\mu\_0}{4\pi} \int\_{V\_\varepsilon} \frac{\mathbf{J}(\mathbf{r}',t)}{|\mathbf{r}-\mathbf{r}'|} dV' + \mathbf{A}\_\varepsilon(\mathbf{r},t) \tag{3}$$

and

$$\mathbf{A}\_{\boldsymbol{\epsilon}}(\mathbf{r},t) = \frac{\mu\_0}{4\pi} \int\_{\mathbb{R}^3 - V\_{\boldsymbol{\epsilon}}} \frac{\mathbf{J}\_{\boldsymbol{\epsilon}}(\mathbf{r}',t)}{|\mathbf{r} - \mathbf{r}'|} dV'. \tag{4}$$

By imposing the Ohm law in weak form, we have

$$\int\_{V\_{\varepsilon}} \mathbf{W} \cdot \left\{ \sigma^{-1} \mathbf{J} + \frac{\partial}{\partial t} \left[ \frac{\mu\_0}{4\pi} \int\_{V\_{\varepsilon}} \frac{\mathbf{J}(\mathbf{r}', \ t)}{|\mathbf{r} - \mathbf{r}'|} dV' \right] + \frac{\partial \mathbf{A}\_{\varepsilon}}{\partial t} \right\} dV = \mathbf{0} \tag{5}$$

for any **W**∈*S* and with **J**∈ *S*, *S* being the set of solenoidal vector functions with continuous normal component across *VC* and zero normal component on the boundary of *VC*.

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 **J** as **J r**ð Þ¼ , *<sup>t</sup>* <sup>P</sup> *<sup>k</sup>Ik*ð Þ*t* ∇ � **N***k*ð Þ**r** , where **N***<sup>k</sup>* 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:

$$\begin{cases} \mathbf{L}\frac{d\mathbf{I}}{dt} + \mathbf{R}\mathbf{I} = -\frac{d\mathbf{V}\_0}{dt}, & t \ge 0 \\\ \mathbf{I}(0) = \mathbf{i}\_0 \end{cases},\tag{6}$$

where **<sup>I</sup>**ðÞ¼ *<sup>t</sup>* ½ � *Ii*ð Þ*<sup>t</sup>* ,**<sup>R</sup>** <sup>¼</sup> *Rij* � �,**<sup>L</sup>** <sup>¼</sup> *Lij* � �,**V**0ðÞ¼ *<sup>t</sup>* ½ � *<sup>V</sup>*0,*<sup>i</sup>*ð Þ*<sup>t</sup>* and

$$R\_{\vec{\eta}} = \int\_{V\_{\varepsilon}} \nabla \times \mathbf{N}\_{\mathrm{i}}(\mathbf{r}) \cdot \sigma^{-1} \nabla \times \mathbf{N}\_{\mathrm{j}}(\mathbf{r}) dV,\tag{7}$$

$$L\_{ij} = \frac{\mu\_0}{4\pi} \int\_{V\_\varepsilon} \int \frac{\nabla \times \mathbf{N}\_i(\mathbf{r}) \cdot \nabla \times \mathbf{N}\_j(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} dV dV',\tag{8}$$

$$V\_{0,k}(t) = \int\_{V\_\varepsilon} \nabla \times \mathbf{N}\_k(\mathbf{r}) \cdot \mathbf{A}\_\varepsilon(\mathbf{r}, t)dV. \tag{9}$$

The quantity **i**<sup>0</sup> in (6) denotes the prescribed initial condition. The vector potential **A***s*ð Þ **r**, *t* is produced by a known source current density **J***s*ð Þ **r**, *t* , and it can be computed by means of the Biot-Savart law. For the sake of simplicity, we assume that the source current density **J***<sup>s</sup>* is due to the current *is*ð Þ*t* flowing in a given coil.

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

$$\begin{cases} \mathbf{AI}\_n = \mathbf{LI}\_{n-1} + \mathbf{V}\_{0,n} - \mathbf{V}\_{0,n-1} \\ \mathbf{I}\_0 = \mathbf{i}\_0 \end{cases},\tag{10}$$

where **I***<sup>n</sup>* ¼ **I**ð Þ *n*Δ*t* and **V***<sup>n</sup>* ¼ **V**ð Þ *n*Δ*t* , and matrix **A** is expressed as

$$\mathbf{A} = \mathbf{R}\Delta\mathbf{t} + \mathbf{L}.\tag{11}$$

The quantity Δ*t* is the time integration step. Matrix **R** is sparse and local, whereas **L** is fully populated because it arises from the discretization of an integral operator. Note that **R**, **L**, and **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 **P** built upon the H<sup>2</sup> -matrix representation of matrix **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 **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 **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{A}\mathbf{I}\_n = \mathbf{P}^{-1}(\mathbf{L}\mathbf{I}\_{n-1} + \mathbf{V}\_{0,n} - \mathbf{V}\_{0,n-1}),\tag{12}$$

where the preconditioned stiffness matrix **P**�**<sup>1</sup> 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 **A**�<sup>1</sup> . A critical issue in this context is the dependence of matrix **A** from Δ*t*. In fast transient analysis, Δ*t* has to be small enough to model properly the underlying dynamics. On the other hand, in slow transients Δ*t* has to be large enough to cover the full-time interval of interest at acceptable computational cost. The dependence of **A** on Δ*t* yields the dependence of **P** on Δ*t*. According to (8), for small Δ*t* the preconditioner should be tailored mainly on **L**, whereas for large Δ*t* the preconditioner should be tailored mostly on **R**. Here and in the following, "small" and "large" Δ*t* should be intended as compared with the slowest eigenmodes of the dynamical matrix **L**�<sup>1</sup> **R**. In intermediate cases (Δ*t* neither too small nor too large), the preconditioner should depend on both **R** *An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

**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*© *code (red), under the hypothesis of axisymmetry.*

and **L**. Since the use of an implicit time stepping does not impose specific restrictions on Δ*t*, the appropriate value of Δ*t* depends on the underlying dynamics induced by the source, i.e., on the rate of change of **V**0ð Þ*t* defined in (9) and appearing in (6).

A canonical waveform for *is*ð Þ*t* , which can be used as the elementary building block to approximate more complex waveforms<sup>1</sup> , is given by a linear ramp *is*ðÞ¼ *t αt*, for *t*≥0, being *α* the rate of change for the electrical current. Consequently, **V**0ðÞ¼ *t* **W***t* for *t*≥0 is also linear. The source term **V**0,*<sup>n</sup>* � **V**0,*n*�<sup>1</sup> is equal to **W**Δ*t*, which is constant for any *n*. Let us consider a transient starting from a vanishing initial condition **i**<sup>0</sup> ¼ **0** in (6). As long as the current at the previous step is negligible, Eq. (10) reduces to **AI***<sup>n</sup>* ¼ **W**. Therefore, the matrix to be inverted is **A**, and the preconditioner **P** has to be "tuned" on such matrix. On the other hand, when the stationary state is achieved, the solution **I***<sup>n</sup>* is constant, and Eq. (10) reduces to ð Þ Δ*t***R** þ **L I**<sup>∞</sup> ¼ **LI**<sup>∞</sup> þ **W**, that is, Δ*t***RI**<sup>∞</sup> ¼ **W**. In this case, the matrix to be inverted is **R**, and the preconditioner needs to be "tuned" on **R**. This explains the reason why the preconditioner based on matrix **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 *is*ðÞ¼ *<sup>t</sup> <sup>α</sup>t*, for *<sup>t</sup>*≥0, where *<sup>α</sup>* <sup>¼</sup> 107*Hz*. This electrical current flows in a circular ring (*R* ¼ 10*cm*), centered on the axis of a flat conducting disk (radius *R*<sup>0</sup> ¼ 1*m*, thickness 3*mm*). Validation is carried out against the Femlab© code [Version 3.51], as shown in **Figure 1** (right). **Figure 1**(left) shows the induced current distribution at *t* ¼ 3*:*4 *ms* together with the finite element mesh used in the computation and the ohmic losses in the conductor. **Figure 2** shows the characteristics

**Figure 2.** *Number of iterations of GMRES using R as a preconditioner during the transient.*

<sup>1</sup> The linear ramp can be used to build a piecewise linear approximation of a continuous function.

of the iterative inversion using **P** ¼ **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 ð Þ **P** ¼ **R** is not effective at the early stage of the transient phase, when the inductive effects accounted by matrix **L** are not negligible.

## **4.** H**2-matrix representation**

The H<sup>2</sup> -matrix representation is widely used for storing efficiently a dense matrix arising from the discretization of a boundary integral equation [15]. In H<sup>2</sup> -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

**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.*

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

#### **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.*

contained within a parallelepiped, which is parallel to the coordinate planes (see **Figure 5**). The key parameters related to a cluster *Nk* are: (i) the number of DoFs contained in the cluster, (ii) the parallelepiped *Pk* containing the DoFs of *Nk*, (iii) the diameter *dk* ¼ j j *Pk* and the center *ck* of *Pk*.

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 *Nk* and *Nj* are named far if the ratio

$$\xi\_{kj} \hat{=} \frac{\left| \mathbf{c}\_k - \mathbf{c}\_j \right|}{\max \left( d\_k, \, d\_j \right)},\tag{13}$$

is greater than a suitable threshold *ξTHR* >1*:*

The mesh partitioning into clusters of nodes is represented by means of the socalled 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 **L** and, hence, of **A** into small blocks, yielding naturally a decomposition **A** ¼ **A***near* þ **A***far*.

Interactions of near clusters lead to high or full-rank blocks (the near part **A***near* of **A**) that need to be explicitely computed and stored, whereas far clusters interactions give rise to low-rank blocks (the far part **A***far* of **A**). In the H<sup>2</sup> -matrix representation, the far interactions between, say, the *i*-th and *j*-th clusters, are approximated by a lowrank matrix factorization of the form

$$\mathbf{A}\_{\vec{\text{ij}}} \approx \mathbf{V}\_i \cdot \mathbf{S}\_{\vec{\text{ij}}} \cdot \mathbf{V}\_j^{\text{T}}.\tag{14}$$

**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.*

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 **V***<sup>i</sup>* represent the expansion basis for the DoFs associated to the *i*-th cluster, and **S***ij* is the coupling matrix of dimensions *ki* � *kj* where *ki* and *kj* are the column dimensions of **V***<sup>i</sup>* and **V***j*, respectively. It is worth noting that for far interactions, integers *ki* and *kj* are expected to be much smaller than the number of rows of matrices **V***<sup>i</sup>* and **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 (**V***<sup>i</sup>* for any *i*, and **S***ij* for any *i* and *j*). Indeed, linear memory complexity can be achieved by means of the H<sup>2</sup> -matrix representation [12].

Hereafter, the following notation will be used: we denote by *Np* the number of Lagrange points per direction used for the degenerate kernel approximation, and by *eij* the relative error of the low-rank representation of a matrix block **A***ij* associated to the interactions of the *i*-th and *j*-th far clusters:

$$e\_{ij} = \frac{||\mathbf{A}\_{ij} - \mathbf{V}\_i \mathbf{S}\_{ij} \mathbf{V}\_j^{\mathrm{T}}||}{||\mathbf{A}\_{ij}||},\tag{15}$$

where ∥ � ∥ is the conventional 2-norm for matrices.

## **5. The** H**2-matrix factorization**

To develop an effective preconditioner **P** for solving Eq. (12), first we set **P** equal to **A**H<sup>2</sup> , the H<sup>2</sup> -matrix representation of **A**. Then, by exploiting the hierarchical structure of **A**H<sup>2</sup> , we build an efficient LU factorization of **P**. A core operation of the factorization of **P** is the computation of the Shur complement for the block matrix **A**H<sup>2</sup> . In general, this computation is very expensive. However, the H<sup>2</sup> -approximation introduces many off-diagonal blocks, which are vanishing upon a suitable orthonormal transformation of **A**H<sup>2</sup> . Hence, the overall factorization costs of **P** can be significantly reduced to Oð Þ *N* � log ð Þ *N* , *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

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

reference to the leaf level, the last level in the BCT, the original stiffness matrix **P** is expressed as:

$$\mathbf{P} = \mathbf{L\_{o}} \mathbf{C\_{1}} \mathbf{U\_{p}} \tag{16}$$

where **C**<sup>1</sup> ¼ *diag* **I***<sup>d</sup>*<sup>1</sup> , **I***<sup>d</sup>*<sup>2</sup> , … , **I***<sup>d</sup>*#ð Þ *leaf clusters* , **C** and **Lo** and **Up** can be easily inverted because they are lower and upper triangular matrices, respectively. In (16), **I***di* is an identity matrix of proper size, and **C** is a full matrix of dimension much smaller than the original matrix **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 **C**.

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

Specifically, matrix **P** is factorized in terms of a partial LU decomposition [16]. For doing this operation efficiently, using some properties of H<sup>2</sup> -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 **A***curr* ¼ **A**H<sup>2</sup> .

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

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

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

5. Compute the partial LU factorization of **A***curr* as.

**A***curr* ¼ **L***<sup>i</sup>* **I 0 <sup>0</sup> <sup>A</sup>**<sup>~</sup> *curr* **<sup>U</sup>***i*, where **<sup>L</sup>***<sup>i</sup>* and **<sup>U</sup>***<sup>i</sup>* are lower and upper triangular matrices, respectively.

$$\begin{array}{cc} \mathbf{\tilde{O}}. & \mathbf{Set} \, \mathbf{A}\_{curr} = \begin{bmatrix} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{\tilde{A}}\_{curr} \end{bmatrix}. \end{array}$$

7. **end for.**

8.**end for.**

At step 5, matrix **Q***<sup>i</sup>* encompasses the basis related to the *i*-th cluster as follows:

**<sup>Q</sup>***<sup>i</sup>* <sup>¼</sup> *diag* **<sup>I</sup>***<sup>d</sup>*<sup>1</sup> , **<sup>I</sup>***<sup>d</sup>*<sup>2</sup> , … , **<sup>V</sup>**<sup>⊥</sup> *<sup>i</sup>* , **V***<sup>i</sup>* , … , **<sup>I</sup>**#ð Þ *leaf clusters :* (17)

The orthonormal transformation at step 5 is computationally cheap to carry out because of the simplified structure of matrix **Q***i*, which involves only the *i*-th diagonal block of matrix **<sup>A</sup>***curr*, of dimension *di* � *di*. In the definition of **<sup>Q</sup>***<sup>i</sup>* in Eq. (17), **<sup>V</sup>**<sup>⊥</sup> *<sup>i</sup>* is the orthogonal complement of the cluster basis **V***i*. This update of **A***curr*, which essentially corresponds to a change of basis, is very effective on the far blocks to the *i*-th diagonal block of **A***curr*. Thanks to (14), it is easy to verify that the first *di* � *k*<sup>ℓ</sup> rows of these blocks are zeroed out, where *k*<sup>ℓ</sup> is the rank used for the approximate factorization at the ℓ-th level of the BCT. Since most of the blocks in the original matrix **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 **A**~ *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 **A***curr* writes as

$$\mathbf{A}\_{curr} = \begin{bmatrix} \mathbf{A}\_{i'i'} & \mathbf{A}\_{i'p} \\ \mathbf{A}\_{pi'} & \mathbf{A}\_{pp} \end{bmatrix} = \underbrace{\begin{bmatrix} \mathbf{I}\_{r'i'} & \mathbf{0} \\ \mathbf{A}\_{r'i'}\mathbf{U}\_{i'}^{-1} & \mathbf{I} \end{bmatrix} \begin{bmatrix} \mathbf{I} & \mathbf{O} \\ \mathbf{O} & \tilde{\mathbf{A}}\_{curr} \end{bmatrix}}\_{\mathbf{L}\_i} \underbrace{\begin{bmatrix} \mathbf{U}\_{i'i'} & \mathbf{I}\_{p i'}^{-1}\mathbf{A}\_{p'p} \\ \mathbf{A}\_{p'i} & \mathbf{I} \end{bmatrix}}\_{\mathbf{U}\_i} \tag{18}$$

where **A***<sup>i</sup>* 0 *i* <sup>0</sup> ¼ **L***<sup>i</sup>* 0 *i* <sup>0</sup>**U***<sup>i</sup>* 0 *i* <sup>0</sup>, *i* <sup>0</sup> is equal to the number of nodes in the *i*-th clusters minus the approximation rank at the current level, and

$$
\tilde{\mathbf{A}}\_{curr} = \mathbf{A}\_{pp} - \mathbf{A}\_{pi'} \mathbf{U}\_{i'i'}^{-1} \mathbf{L}\_{i'i'}^{-1} \mathbf{A}\_{i'p}.\tag{19}
$$

Clearly, matrix **A**~ *curr* still has a block structure, which is inherited from the clusterto-cluster interactions. Contributions to the update term **A***pi*0**U**�<sup>1</sup> *i* 0 *i* <sup>0</sup> **L**�<sup>1</sup> *i* 0 *i* <sup>0</sup> **A***<sup>i</sup>* 0 *<sup>p</sup>* 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 **V***i*, thus avoiding to re-compute the whole basis from scratch. This approximate computation of the Schur complement is essential to achieve Oð Þ *N* � log ð Þ *N* complexity. For further details, see [16].

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


### **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**.

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

**Figure 7.** *The analyzed case study. Left: mesh of ITER magnet system. Right: A* 45<sup>∘</sup> *sector.*

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).

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 Δ*t* is set equal to 1*ms* and 0*:*1*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 **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 **A**H<sup>2</sup> 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 **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 *Np* and the leaf size *LS*. We remind that *LS* determines the size of matrix **A***near* arising from near cluster interactions, whereas *Np* affects the number of columns of **V***i*, the size of **S***ij* and, consequently, the accuracy of **A***far*. It is worth noting that **A***far* in the H<sup>2</sup> -matrix representation does not contribute significantly to the overall required memory: in all numerical experiments of this work, the memory required by **A***far* is less than 10% of the storage required by **A***near*.

From (11), we have that the preconditioner memory footprint depends only on the storage of matrices **L**0, **U***<sup>p</sup>* and the dense but smaller block **C** of matrix **C**<sup>1</sup> appearing in (16). **Figure 8** shows the memory required by the preconditioner (*MLU*) as a function of *LS* and *Np*. The result is normalized with respect to the memory required to store the full matrix **A**, which is about 169 *GB*.

The memory footprint required by matrix **C**<sup>1</sup> decreases with *LS*. Indeed, an increase of *LS* enlarges the size of the block-to-block interactions, resulting in a smaller reduced matrix **C**<sup>1</sup> because of a more effective compression of the block-toblock interactions. On the other hand, the sizes of **L**<sup>0</sup> and **U***<sup>p</sup>* depend on: (i) the size of the near part, strictly depending on *LS*, and (ii) the ability of the H<sup>2</sup> -matrix representation to compress the admissible block-to-block far-field interactions.

If we reduce *LS*, the memory required by the near interactions matrix **A***near* decreases, but the memory required by the far-field block-to-block interactions matrix increases. Indeed, we observe that at *LS* ¼ 500 the required memory depends on *Np*, whereas for *LS* ¼ 2000 the required memory is almost independent on *Np*. Eventually, we have an intermediate dependence for *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 **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 **R** plus a contribution from matrix **L** on the same pattern of **R**. This preconditioner failed to converge at increasing mesh size.

**Figure 8.** *Memory occupation as LS and Np varies.*

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

Summing up, finding a proper preconditioner for this class of large-scale problems is challenging.

**Figure 9** (top Δ*t* ¼ 1*ms*) and (bottom Δ*t* ¼ 0*:*1*ms*) shows the preconditioned relative residual norm *ER* of the GMRES solver versus the number of iterations for different *LS* and *Np*.

The numerical results clearly show that:

• The H<sup>2</sup> -LU generally performs much better than the preconditioner based on **R** preconditioner.

**Figure 9.**

*Performance of* H<sup>2</sup> *-LU preconditioner, for different values of the leaf size (LS) and number of Lagrange points (Np). Top:* Δ*t* ¼ 1 *ms. bottom:* Δ*t* ¼ 0*:*1 *ms. performances from the reference preconditioner based on R are also shown.*

• The efficiency of the preconditioner strongly depends on the choice of the parameters *Np* and *LS* that control the accuracy of **A***far*.

We highlight that the accuracy of **A***far* is rather low for *Np* ¼ 1 and that the preconditioner based on **R** performs better in this situation. In the other cases, the proposed H<sup>2</sup> -LU preconditioner outperforms the reference one based on **R**, proving that it incorporates a good mechanism to balance contributions from matrices **L** and **R**Δ*t* in **A**. This explains why at the smaller Δ*t*, when **L** gives the dominant contribution to **A**, the number of iterations for the **R** based preconditioner increases. On the other hand, when using the H<sup>2</sup> -LU preconditioner, which accounts for both **R**Δ*t* and the dominant part of **L**, the number of iterations remains almost unchanged with Δ*t*.

The accuracy of **A***far* is a critical point for the effectiveness of the proposed preconditioner. To asses this, we define two quantities, denoted as *EA* and *ESD*, that are the mean and the standard deviation for the relative error *eij* defined in (15) arising from all the admissible block-to-block pairs ð Þ *i*, *j* . As expected, **Figure 10** shows that the error decreases as *Np* increases. Note that for *Np* ¼ 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 *LS*, *Np*, and Δ*t*. To quantify the error, we introduce the metric *ES* defined as

$$E\_S = \frac{||I\_{app} - I||}{||I||},\tag{20}$$

where *I* is computed using a direct method, whereas *Iapp* is computed with the proposed preconditioned iterative solver. In (19), k k� is the 2-norm, and both *I* and *Iapp* are evaluated at a prescribed time.

**Figure 10.** *Mean and standard deviation of the relative error in Afar.*

**Figure 11.** *Solution error Es, as LS, Np, and* Δ*t vary.*

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

**Figure 11** confirms that the accuracy of the solution increases with *Np*. Only for *LS* ¼ 500 and Δ*t* ¼ 1*ms* the solution with *Np* ¼ 3 is a little bit less accurate than the one with *Np* ¼ 2. It is worth noting that as the accuracy of **A***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 **L**.

### **7. Conclusions and future work**

In this paper, a new preconditioner based on the H<sup>2</sup> -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 H<sup>2</sup> -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.

### **Acknowledgements**

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.

## **Author details**

Salvatore Ventre<sup>1</sup> , Bruno Carpentieri<sup>2</sup> \*, Gaspare Giovinco<sup>3</sup> , Antonello Tamburrino1,4, Fabio Villone<sup>5</sup> and Guglielmo Rubinacci<sup>5</sup>

1 Department of Electrical Engineering and Information, University of Cassino and Southern Lazio, Cassino, Italy

2 Faculty of Computer Science, Free University of Bozen-Bolzano, Bolzano, Italy

3 Department of Civil and Mechanical Engineering, University of Cassino and Southern Lazio, Cassino, Italy

4 Department of Electrical and Computer Engineering, Michigan State University, East Lansing, USA

5 Department of Electrical Engineering and Information Technology, University of Naples Federico II, Napoli, Italy

\*Address all correspondence to: bruno.carpentieri@unibz.it

© 2022 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*An Effective* H2 *-LU Preconditioner for Iterative Solution of MQS Integral-Based Formulation … DOI: http://dx.doi.org/10.5772/intechopen.108106*

## **References**

[1] Albanese R, Rubinacci G. Finite element methods for the solution of 3D eddy current problems. Advances in Imaging and Electron Physics. 1997;**102**: 1-86

[2] Rubinacci G, Tamburrino A, Villone F. Circuits/fields coupling and multiply connected domains in integral formulations. IEEE Transactions on Magnetics. 2002;**38**(2):581-584

[3] Rubinacci G, Tamburrino A. Automatic treatment of multiply connected regions in integral formulations. IEEE Transactions on Magnetics. 2010;**46**(8):2791-2794

[4] Carpentieri B. Algebraic preconditioners for the fast multipole method in electromagnetic scattering analysis from large structures: Trends and problems. Electronic Journal of Boundary Element. 2009;**7**(1):13-49

[5] Carpentieri B. Fast iterative solution methods in electromagnetic scattering. Progress in Electromagnetic Research. 2007;**79**:151-178

[6] Carpentieri B. A matrix-free two-grid preconditioner for boundary integral equations in electromagnetism. Computing. 2006;**77**(3):275-296

[7] Carpentieri B, Bollhöfer M. Symmetric inverse-based multilevel ILU preconditioning for solving dense complex non-hermitian systems in electromagnetics. Progress In Electromagnetics Research. 2012;**128**: 55-74

[8] Carpentieri B, Duff I, Giraud L, Sylvand G. Combining fast multipole techniques and an approximate inverse preconditioner for large electromagnetism calculations. SIAM

Journal on Scientific Computing. 2005; **27**(3):774-792

[9] Albanese R, Rubinacci G. Solution of three dimensional eddy current problems by integral and differential methods. IEEE Transactions on Magnetics. 1988;**24**(1):98-101

[10] Rubinacci G, Ventre S, Villone F, Liu Y. A fast technique applied to the analysis of resistive wall modes with 3D conducting structures. Journal of Computational Physics. 2009;**228**(5): 1562-1572

[11] Hackbusch W. A sparse matrix arithmetic based on ℋ-matrices. Part i: Introduction to ℋ-matrices. Computing. 1999;**62**(2):89-108

[12] Börm S, Grasedyck L, Hackbusch W. Hierarchical Matrices. Leipzig: Max Planck Institute for Mathematics in the Sciences; 2003

[13] Voltolina D, Bettini P, Alotto P, Moro F, Torchio R. High-performance PEEC analysis of electromagnetic scatterers. IEEE Transactions on Magnetics. 2019;**55**(6):1-4

[14] Bebendorf M. Hierarchical LU decomposition-based preconditioners for BEM. Computing. 2005;**74**(3): 225-247

[15] Börm S. Efficient Numerical Methods for Non-Local Operators: ℋ2- Matrix Compression, Algorithms and Analysis. Vol. 14. Zürich, Switzerland: European Mathematical Society; 2010

[16] Ma M, Jiao D. Accuracy directly controlled fast direct solution of general ℋ2-matrices and its application to solving Electrodynamic volume integral equations. IEEE Transactions on

Microwave Theory and Techniques. 2017;**66**(1):35-48

[17] Albanese R, Rubinacci G. Integral formulation for 3D eddy-current computation using edge elements. IEE Proceedings A-Physical Science, Measurement and Instrumentation, Management and Education-Reviews. 1988;**135**(7):457-462

[18] Atkinson K, Han W, Stewart D. Numerical Solution of Ordinary Differential Equations. Vol. 108. Hoboken, New Jersey, U.S.: John Wiley & Sons; 2011

[19] Golub G, Loan CV. Matrix Computations: Johns Hopkins Studies in the Mathematical Sciences. third ed. Baltimore, MD, USA: The Johns Hopkins University Press; 1996

[20] Saad Y. Iterative Methods for Sparse Linear Systems. 2nd ed. Philadelphia: SIAM; 2003

[21] Ventre S, Carpentieri B, Giovinco G, Tamburrino A, Rubinacci G, Villone F. An ℋ2-LU preconditioner for MQS model based on integral formulation. In: Poster Presentation at the 22nd International Conference on the Computation of Electromagnetic Fields, COMPUMAG 2019, July 15–19. Paris, France; 2019

[22] Testoni P, Forte R, Portone A, Rubinacci G, Ventre S. Damping effect on the ITER vacuum vessel displacements during slow downward locked and rotating asymmetric vertical displacement events. Fusion Engineering and Design. 2018;**136**:265-269

[23] Ramogida G, Maddaluno G, Villone F, et al. First disruption studies and simulations in view of the development of the DEMO physics basis. Fusion Engineering and Design, Open Access. 2015;**96–97**:348-352

[24] Ariola M, Pironti A. Magnetic Control of Tokamak Plasmas. Vol. 187. Berlin/Heidelberg, Germany: Springer; 2008

[25] Hender T, Wesley J, Bialek J, et al. MHD stability, operational limits and disruptions. Nuclear Fusion. 2007;**47**(6): S128

## **Chapter 3**
