An efficient two-level preconditioner based on lifting for FEM-BEM equations

An efficient two-level preconditioner based on lifting for FEM-BEM equations

Fabio Henrique Pereira^I; Marcio Matias Afonso^II; João Antônio de Vasconcelos^III; Silvio Ikuyo Nabeta^IV

^IIndustrial Engineering Post Graduation Program – UNINOVE University-São Paulo, Brazil

^IIElectrical Engineering Department - Federal Technology Education Center - CEFET-Minas Gerais, Brazil

^IIIElectrical Engineering Department – UFMG - Minas Gerais, Brazil

^IVElectrical Machine and Drives Lab – São Paulo University, GMAcq/USP – São Paulo, Brazil

ABSTRACT

The system resulting from the coupled Finite Element Method and Boundary Element Method formulations inherits all characteristics of both finite element and boundary element equation system, i. e., the system is partially sparse and symmetric and partially full and nonsymmetric. Consequently, to solve the resulting coupled equation system is not a trivial task. This paper proposes a new efficient lifting-based two level preconditioner for the coupled global system. The proposed approach is applied to solve the coupled systems resulting from the electromagnetic scattering problem and its performance is evaluated based on the number of iterations and the computational time. Traditional methods based on incomplete and complete LU decompositions are used for comparison.

Index Terms: FEM-BEM Equations, Lifting Technique, Preconditioner.

I- INTRODUCTION

The Finite Element and the Boundary Element Methods can be considered complementary to each other in many cases. It is particularly true to problems where the domain considered is formed by two (or more) regions where one of them is characterized by non-homogeneous material while the other is free space. This kind of problem is difficult to be solved by any individual method described above. Therefore, a combination of the two previous techniques comes up as a powerful way to solve it [1].

The system resulting from the coupled Finite Element Method and Boundary Element Method (FEM–BEM) formulations inherits all characteristics of both finite element and boundary element equation system, i. e., the system is partially sparse and symmetric and partially full and nonsymmetric [1], [2]. Consequently, the efficient solution of the FEM-BEM system has been object of many researches [3], [4]. In fact, most of the proposed approaches are based on some adaptations. In the algebraic multigrid context, for example, one should consider a fast matrix-by-vector multiplication for dense matrices and to use some matrix approximation technique for the BEM dense matrix. In particular, common smoother types do not work properly in the case of the single layer potential [5]. In face of that, this work proposes a new efficient two level preconditioner for the global FEM-BEM system based on the lifting technique [13]. This proposed approach is applied to some real-life FEM-BEM coupled systems and its performance is evaluated based on the number of iterations and the computational time. Traditional methods based on incomplete and complete LU decompositions are used for comparison.

II. PROBLEM FORMULATION

The scattering problem was chosen to show the performance of the proposed approach. This kind of problem occurs when an object is struck by an electromagnetic wave. To solve this problem by hybrid technique an artificial surface must be chosen, which allows dividing scattering problems in two regions. The first, , is the free space with permeability µ₀ and permittivity ε₀. The second region Ω may consist in general of non-homogeneous material with permeability µ(x, y) and permittivity ε(x, y). Once the two regions are known, the electromagnetic fields in each one can be formulated. These formulations for both interior and exterior fields can be then coupled at the boundary surface through continuity conditions [1],[2],[6]. In this work, the formulation for both TMz and TEz polarization are provided, the e^jwttime convention is assumed and, the incident field is given by

A. Finite element formulation

The Helmholtz differential equation describes the field behavior in region Ω. In general form, the Helmholtz equation could be written as [1]

where, represents the wave number. Also, for electric field polarization , and u=E_z while for magnetic field polarization , and u = H_z.

The strong problem formulation for this case can be expressed as it follows: given α₁, α₂ and k₀, find u such that [1]:

In above equations , and are exterior and Dirichlet boundary, respectively. Also, ψ is the normal derivative and is the Dirichlet potential.

The weak formulation given by (6) is well known and its development will be omitted here [2].

in which u is an approximation for the field and w is the weight function. The finite element formulation given by (6) is computed in the target domain and it is well suited to deal with complex geometries and anisotropic or isotropic materials [2].

Applying Galerkin's procedure with (6) and dividing the domain in a grid of triangular elements, it is possible to determine the general set of equations that can be written in matrix form as [2],

where [K] is a n x n square matrix and [C] is a n x m rectangular matrix, with m and n representing the total number of the nodes on the boundary and in the interior domain, respectively. Also, {d} is a column vector for approximated field arguments u and {k} is a column vector for the normal derivative field arguments.

B. Boundary element formulation

In free space, , the fields also are formulated by Helmholtz equation. The general formulation for both electric and magnetic fields is:

in which, for electric polarization u=E_z, and for magnetic polarization u = H_z, . Where, Z₀ is the intrinsic impedance of the free space, J_z represents all currents that flow along z-axis and M_i is the impressed magnetic source.

To formulate the fields in , one introduces the free space Green's function G₀, which satisfies the differential equation

and the Sommerfeld radiation condition [2]. In Eq. (9), represents the Dirac delta function. The solution for this equation is [2]

The boundary fields are derived from the product of Eq. (8) by G₀, the application of second scalar Green's theorem and Dirac properties [1],[2].

In Eq. (11), the primed sign indicates operation under integration point, the factor 0.5 is the solid angle saw by the observer on the boundary and is the incident field given by (1).

Then, writing Eq. (11) for all m nodes on the boundary, it produces [1],[2]

in which, [H] and [Q] are m x m matrices obtained from the discretization and {b} is a m x 1 incident field column vector.

C. Fem-Bem formulations

A typical finite element – boundary element system consists of two sets of equations, where one of them is obtained from FEM (7) and other obtained on the boundary from BEM (12). As the unknowns in both equations are the same, the two sets of equations could be coupled and expressed in the in following form:

This equation system inherits all characteristics of both finite element and boundary element methods. For this reason, the equation system given by (13) is partially sparse and symmetric and partially full and nonsymmetric [1]. There are many different ways to write FEM–BEM equation, however the form present here is the easiest one [2]. The system is composed by N x N algebraic equations, N = (m + n), and solves both finite element and boundary element equation systems simultaneously. Although (13) is considered the best form to understand the coupled FEM–BEM method it is not the better way to solve it because the symmetry of the FEM subsystem is not exploited [2]. To write this system more efficiently see [2], [7]. Especial attention is required to evaluate [H] and [Q] because there are singularities when both source and observation points coincide [1], [2].

III-Lifting Based Two-Level Preconditioner

The approach proposed here is based on a projection method. In this process an approximation for the solution of the linear system (13) is extracted from a subspace of which is the search subspace. If p is the dimension of , then the residual vector b- Ax is constrained to be orthogonal to the same subspace . Mathematically, denoting by V a N x p matrix whose column-vectors form a basis of , with an initial guess , the approximate solution is written as

in which y is obtained solving the following system of equations from the orthogonality condition [8]

in which the initial residual vector is r₀ = b - Ax₀.

Then, if the p x p matrix V^T AV is nonsingular, the expression for the approximate solution becomes,

In the lifting based approach, the matrix V^T AV is obtained by the decomposition of the matrix A induced by the discrete lifting wavelet transform, with the wavelet of Haar [9].

The multiresolution analysis in discrete lifting wavelet transform decomposes the original grid (space) Ω ∪ in two subspaces V and W such that

and

whereand ψ are the scaling and wavelets function, respectively, V is the low frequency and W the high frequency subspaces. The direct sum V ⊕ W in (19) means the sum of subspaces in which V and W have only the zero vector in common [10].

Therefore, this approach can be viewed as a two-level multigrid algorithm in which V^T is the restriction operator that represents a wavelet transform in equation form (lifting), as proposed in [11]. There are basically three forms for representing a wavelet transform: equation form (lifting), filter form (filter bank) and matrix form. However, only the first two are appropriate in the multigrid implementation. In filter form, the restriction operator V^T can be defined as (20), for an original system with N unknowns,

in which are the low-pass filters coefficients associated to the haar wavelet.

In this work, however, the application of the lifting restriction operator V^T is algebraically defined using the equation form of the wavelet transform as stated in (21), in which s represents an approximation in the coarse grid of the original vector S, i. e., s = V^TS.

The main advantages of this technique over the classical wavelet transform are:

a) Smaller memory requirement – the calculations can be performed in-place;

b) Efficiency: reduced number of floating point operations;

c) Parallelism –inherently parallel feature;

d) Easier to understand - not introduced using the Fourier transform;

e) Easier to implement;

f) The inverse transform is easier to perform – it has exactly the same complexity as the forward transform;

g) Transforms signals with a finite length (without extension of the signal).

More details about these advantages as well as other important structural advantages of the lifting can be found in [12], [13].

The computation of the coarse matrix V^T AV is done applying the procedure (21) in the rows and columns of the original matrix A. If the number of rows and columns in matrix A is not even, a zero padding operation is accomplished adding one sample with value zero at the end of the rows and columns in order to achieve the desired length [12]. So, the resulting V^T AV matrix will be formed by the detail coefficients obtained from the discrete lifting wavelet transform and will represent an approximation in the coarse grid of the matrix A which is defined as (13).

The procedure defined in (21) represents the Haar transform in the lifting scheme. Haar is the simplest possible wavelet and it is associated to shorter filters [12]. There are many other wavelet functions that can present better approximation properties in many cases, but the choice for short filters is fundamental to control the fill-in in matrix V^TAV and keep its sparsity, as already discussed in [11].

This two-level preconditioner is illustrated in Fig. 2 and it can be represented by the following algorithm:

As the matrix A is partially sparse and the procedure defined in (21) is the Haar transform, the resulting matrix will also be partially sparse and it can be created explicitly even for very large problems. However, if an iterative method was used to solve the coarse grid system in line 3 of the algorithm the matrix V^T AV has not to be explicitly created.

IV-NUMERICAL RESULTS

The Lifting based two-level method (LTL) was applied as a preconditioner for the Stabilized Bi-Conjugate Gradient Method (BiCGStab) in the solution of three global coupled FEM-BEM systems. These systems were obtained for a circular cylinder with 0,3λ diameter, ε_r = 3 and wave length λ = 1. The cylinder is struck by a TM_z plane wave with incident angle of θⁱ= 180º.

The results related to the number of iterations and setup and solver times are presented in tables I, ^II and ^III. The number of iterations corresponds to the number of BiCGStab steps necessary to reduce the Euclidean norm of the residual vector to the order of 10^-4, which was considered enough to reach the desired accuracy.

In Tables I to III, N and NZ stand for the number of rows and the number of nonzero entries in the matrix, respectively.

In the tests, the system with the V^T AV matrix in LTL was solved by incomplete LU solver (LTL-1) and with a direct method based on the complete LU decomposition (LTL-2). So, the setup time presented includes the computation of a explicit representation of the coarse matrix V^T AV. Such representation of the coarse matrix has a sparsity pattern similar to the original matrix (Fig. 3) and it can be created in a relative small computational time, as can be viewed in the setup time column from tables I, ^II and ^III. The Incomplete LU (ILU) preconditioner and the complete LU solver were used for comparison.

The convergence histories for the tests are illustrated in Fig. 4. The convergence history shows the residual norm after the first preconditioner step what explains the difference in the start point of the convergence curves. As the initial guess was the zero vector, the initial residual norm is equal to the norm of the right hand side vector in (13).

CONCLUSIONS AND COMMENTS

The approach proposed in this work revealed to be very efficient and promising to solve the global systems of coupled FEM-BEM equations, mainly for the medium and large problems. The usual complete LU decomposition method is a good choice for small problems, but is very expensive for larger problems, as can be viewed in ^{Table II}.

For large problems, where an efficient preconditioner is required, the lifting two-level method gives the best results. In this case, the LU decomposition can be used to solve efficiently the coarse system which has the half of unknowns (see ^{Table III}). Concerning the computational time the proposed approach was about 80% faster than the classical ILU preconditioner. Also, regarding the convergence rate, the numerical results suggest that the number of iterations necessary to reach a given accuracy grows only logarithmically with the number of unknowns. Further tests have to be done in order to prove this conjecture.

About the coarse system solution (step 3 in the algorithm), the coarse matrix V^TAV should be explicitly assembled if an incomplete or complete LU solver is used in this step. However, this is not a problem since the resulting matrix is partially sparse as can be seen in Fig. 3, and it can be created explicitly even for very large problems. Also, regarding the memory requested it can be seen that the operator complexity of this preconditioner, which is defined as the sum of number of nonzero entries in the coarse matrix and its ILU factors divided by the number of nonzero entries in the matrix A, is something about 80% of operator complexity in the ILU preconditioner. It indicates the total storage space required by the preconditioner matrices and it is generally considered a good indication of the preconditioning operation cost in multigrid approaches.

One alternative is using an iterative solver such as BiCGStab, where it is enough to know the action of the system matrix on a given vector, but the preconditioning of the coarse system is difficult to be achieved in this case.

Finally, it is important to say that the proposed method can be applied in different kind of linear system, complex or real, symmetric or non-symmetric, sparse or dense matrices, and it can be a good alternative in those cases in which the conventional block-diagonal preconditioning for Schur complement based or hierarchical basis solvers do not produce good results.

ACKNOWLEDGEMENT

The authors would like to thank FAPESP (2009/11082-2) and CNPq for the financial support.

received 3 May 2010

revised 5 August 2010

accepted 5 August 2010

• [1] M. M. Afonso, J. A. Vasconcelos, "Two dimensional scattering problems solved by finite element boundary element method", In: 22th CILAMCE - Iberian Latin-American Congress on Computational Methods in Engineering, 2001.
• [2] J. Jin, The finite element method in electromagnetics, Wiley, New York, 1993.
• [3] G. Aiello, S. Alfonzetti, E. Dilettoso, N. Salerno,  "An Iterative Solution to FEM-BEM Algebraic Systems for Open-Boundary Electrostatic Problems", IEEE Trans. on Magnetics, vol. 43, n. 4, pp. 1249 1252, 2007.
• [4] P. Mund and E. P. Stephan, "Adaptive Coupling and Fast Solution of FEM-BEM Equations for Parabolic-Elliptic Interface Problems", Math. Methods in the Appl. Scienc., vol. 20, n. 5, pp. 403-423, 1998.
• [5] D. Pusch, Efficient Algebraic Multigrid Preconditioners for Boundary Element Matrices, Phd thesis. Institute of Computational Mathematics, Johannes Kepler University Linz, Austria, 1997.
• [6] A. C. Balanis, Advanced engineering electromagnetics, Wiley, New York, 1989.
• [7] N. Lu and J. M. Jin, "Application of fast multipole method to finite element boundary integral solution of scattering problems", IEEE Transaction on antennas and propagation, vol. 44, n. 6, pp. 781 786, 1996.
• [8] Y. Saad, Iterative Methods for Sparse Linear Systems, 2^nd ed., SIAM, Philadelphia, 2003.
• [9] T. K. Sarkar, S. P. Magdalena, C.W. Michael, Wavelet Applications in Engineering Electromagnetics, Artech House, Boston, 2002.
• [10]  K. H. Rosen, (Ed.). Handbook of Discrete and Combinatorial Mathematics. Boca Raton, FL: CRC Press, 2000, pp. 357.
• [11]   F. H. Pereira and S. I. Nabeta, "Wavelet-Based Algebraic Multigrid Method Using the Lifting Technique", Journal of Microwaves, Optoelectronics and Electromagnetic Applications, vol. 9, n. 1, 1-9, 2010.
• [12]  A. Jensen, A. la Cour-Harbo, The Discrete Wavelet Transform, Ripples in Mathematics, Springer, Berlin 2001, pp. 21.
• [13]  W. Sweldens. "The lifting scheme: A custom-design construction of biorthogonal wavelets". Appl. Comput. Harmon. Analysis, 3(2):186-200 1996.

Publication Dates

• Publication in this collection
  08 Nov 2011
• Date of issue
  Dec 2010

History

• Received
  03 May 2010
• Reviewed
  05 Aug 2010
• Accepted
  05 Aug 2010