A Multiresolution Domain Decomposition Preconditioner for the MoM Solution of Multiscale Complex Structures

In this communication, we present the combination of a high scalability implementation of the multibranch–multiresolution preconditioner with the domain decomposition method for the electromagnetic analysis of geometrically complex structures with different levels of multiscale features and discretized with a possible nonconformal mesh. Finally, a numerical experiment is shown to illustrate the great efficiency of the proposed approach for the solution of large multiscale objects.

scenario and slowing down convergence.In these cases, fast iterative solvers usually fail to reach an accurate solution in a reasonable time, and the need for preconditioners to speed up iterative convergence becomes evident.Among many variants available, physics-based preconditioners relying on quasi-Helmholtz decompositions [10], [11] or the Calderón identities [12], [13] have proven good success.Alternatively, algebraic preconditioners, such as incomplete lower unitriangular (ILU) upper triangular, sparse approximate inverse (SPAI), or block-Jacobi [14], [15], [16], [17], [18], also improve convergence considerably.In the case of large-scale problems including multiscale features, Schwarz preconditioners based on the domain decomposition method (DDM) [4], [19], [20], [21], [22] stand out for their significant improvement in convergence, also bringing additional advantages in handling such complex problems.These preconditioners can be categorized as algebraic preconditioners, since they estimate an inverse of the matrix system, and also as physicsbased preconditioners, since they allow to separate the physics of the different subsystems that make up the whole problem to adequately solve each one using the method best tailored to their particular features.
Nonetheless, despite the use of outperforming preconditioners, the effective solution of the increasingly complex high-fidelity models that are demanded by the industry today must be undertaken with extreme care.It is well known that the judicious selection of the geometric partition into subdomains plays a fundamental role in the performance of DDM [6], [23].Choosing many small subdomains can lead to fast local problems, but slow convergence of the external iterative algorithm.On the other hand, choosing larger subdomains will improve the effectiveness of the preconditioner and, thus, the external convergence, but at the expense of a complexity shift toward local solvers, which will face medium-to large-sized ill-conditioned problems.A representative example is the case of a supporting structure containing large yet complex subsystems, with details at multiple scales (e.g., the case of large antenna arrays).Another case is the presence of large cavities that, even if they are geometrically simpler, support strong multiple-bounce interactions.In these and other similar cases, these subsystems must be treated as a whole, seeking not to artificially break the strong internal local interactions.A parallel goal is to isolate as far as possible these strong interactions from the external iterative algorithm, which in this way can be focused on finding the global solution without dwelling on the internal details of the subdomains.
As a result of all the above, the incorporation of efficient preconditioners to the local solvers then becomes essential to improve the convergence in the local stage, resulting in a better overall performance of the DDM approach.Among the different preconditioning alternatives, the multiresolution (MR) preconditioner [11], [24], [25] emerges as a strong candidate to be integrated into a domain decomposition scheme to speed up the solution of the local subdomains.This preconditioner introduces a set of multilevel basis functions to discretize the problem while keeping the different scales of variation of the solution [26], [27], improving the conditioning of the original system [28] through a multilevel quasi-Helmholtz decomposition in which the unknown current is separated into its solenoidal and nonsolenoid parts.The MR preconditioner has proven to be a good choice to improve the convergence in multiscale problems [24], [29], especially in the case of electrically medium-sized geometries, and it can be efficiently embedded into the MLFMA-fast Fourier transform (MLFMA-FFT) framework as a multiplicative preconditioner.In addition, the MR preconditioner is the only multilevel quasi-Helmholtz decomposition able to find automatically topological loops in nonconformal meshes [30].
In this communication, we propose to integrate the multiresolution preconditioner into the DDM to derive an efficient multiscale fast solver capable of addressing large-scale problems including medium-sized complex subsystems exhibiting deep multiscale nature.In addition, the multibranch Rao-Wilton-Glisson (MB-RWG) basis and testing functions [30], [31], [32], [33] are applied.A numerical example is presented to demonstrate the versatility and capability of the proposed method for multiscale electrically large problems.Preliminary results with a theoretical example consisting of an array antenna were recently presented in [34], showing the capabilities of the proposed approach when dealing with nonconformal multiscale objects using the electric field integral equation (EFIE).In this work, the formulation is extended to address complex structures exhibiting multiple levels of detail at different scales, and a combined field integral equation (CFIE)-EFIE (combined EFIEs) formulation is considered for solving closed and open surfaces.Then, the performance of the method is examined through the solution of a high-fidelity realistic model of a satellite with several communication systems, thus demonstrating the capabilities of the proposed scheme for solving cutting-edge applications demanded by the industry.Furthermore, the convergence solution is compared with two efficient preconditioners, commonly applied in the context of a DDM approach.
The rest of this communication is organized as follows.The background of the formulation is briefly reviewed in Section II, setting the notation.In Section III, the formulation of the proposed MR-DDM method is discussed in detail.A numerical result analyzing a challenging problem is presented in Section IV.Finally, some concluding remarks are drawn in Section V.

II. BACKGROUND
We consider the electromagnetic scattering of a perfect electric conductor (PEC) body in a homogeneous background.We define J(r ′ ) as the equivalent electric current density on the body surface.The tangential EFIE (T-EFIE) and the normal magnetic field integral equation (N-MFIE) can be obtained by applying the equivalence principle to the total electric and magnetic fields as follows: where η is the intrinsic impedance of the background, n is the unit vector normal to the surface, and E inc and H inc are the incident electric and magnetic fields, respectively.The integro-differential operators L and K are defined as where k is the wavenumber; r and r ′ are the observation and source point, respectively; ∇ ′ denotes the divergence in the source coordinates; and PV denotes the principal value of the integral in (4).
The homogeneous Green's function g(r, r′) is defined as In order to derive a well-tested formulation, we combine (1) and ( 2), obtaining the CFIE as where 0 < α < 1 (α = 0.5 in this work) is a parameter to balance the weight of the EFIE and MFIE equations.
The MoM procedure is applied to (6), expanding the equivalent electric current densities into a sum of known vector basis functions f n as where I n are the unknown complex coefficients.In (7), both RWG and MB-RWG bases are used, allowing nonconformal triangular meshes at any point of the geometry, including the tear contours between adjacent subdomains, which thereby can be independently generated and meshed.Then, applying the Galerkin testing procedure, we obtain a dense matrix system as follows: where with Z EFIE and Z MFIE equal to N × N matrices, containing the coupling between all the basis and testing functions for the EFIE and MFIE formulations, respectively, and α is the ponderation factor, which is selected equal to 0.5 in this work to optimize the tradeoff between conditioning and accuracy.In (8), [I ] is an N -column vector collecting the unknown coefficients I n of the current expansion (7), and [V ] is the N -column excitation vector, closely related to the incident fields originated by the sources.

III. FORMULATION
Let us start with the matrix system of linear equations posed in ( 8) and ( 9).The original problem can be decomposed into a collection of K subdomains, D k , with k = 1, . . ., K , depending on the geometrical features.Then, an additive Schwarz DDM preconditioner [14] can be applied for the solution of the previous matrix system as a left-preconditioner along the solutions of the individual subdomains as follows: where [P] is the DDM block diagonal preconditioner, which can be written as Each diagonal block P k formally represents the impedance matrix of the respective subdomain D k , i.e., P k = Z k .
As detailed in [4], the application of the DDM preconditioner to (10) involves a global matrix-vector product (MVP) representing the global coupling between subdomains and a multiplication by the block diagonal matrix [P] −1 , denoting the local solution of the individual subdomains.These operations are repeated at each iteration of the global (outer) iterative algorithm.Although formally stated with the inverses of the block matrices, the preconditioning local problems of each outer iteration can be solved by the method considered most appropriate in each case.Considering the above, the subdomain problems can be written at each outer iteration as where I k is the local solution and V k is the local right-hand side (RHS), obtained as the global MVP restricted to the D k subdomain.Achieving good convergence with the above preconditioner demands that normal current continuity is satisfied between subdomains in contact.One way to improve this continuity is through the concept of domain enlarging, using the so-called tear-andinterconnect DDM approach [4].With this technique, the local subdomains solved at each outer iteration step are enlarged to incorporate the near-field current flowing across the tearing contours.The enlargement is done by including "flaps" [4] of a quarter to a half wavelength width, which indeed belongs to the adjacent touching subdomains.The reader interested in further application details of DDM can refer to [4] and [21].

A. Local Systems Preconditioning
We focus now on the solution of the local problems stated in (12).To simplify the notation, hereafter, we remove the local indexes denoting subdomains.Hence, without loss of generality, we can rewrite (12) as Notwithstanding, everything that follows in this section applies to the independently solved local system in each subdomain D k .
To expedite (and in some cases enable) the solution of the local matrix system in (13), the MR preconditioner [11] is applied to project the original matrix system into a new function subspace that improves the eigenvalues' distribution.The key point of this effective preconditioner is the application of a quasi-Helmholtz decomposition on the basis of a multilevel scheme.In this work, this decomposition is allowed to start from a nonconformal triangular mesh.The complete set of original RWG and MB-RWG bases is transformed into a set of multilevel bases, generalized basis functions (gF).At each level-l, l = 1, . . ., L, the corresponding gF l values are split to generate the solenoidal and nonsolenoidal basis functions, called latter MR basis, by means of a procedure that automatically includes all the topological solenoidal functions [35].The above scheme is applied recursively down to the quasi-Nyquist cell-size level (level-L), leaving (gF L ) defined at the last level.The obtained multilevel set of basis functions (MR functions of levels l = 1, . . ., L − 1 and gF L ) spans the same space of the original bases [11].
The MoM system matrix can be expressed in the new space of functions as where the matrix [T ] is the change-of-basis matrix, collecting the coefficients of the new basis functions written as a linear combination of the initial RWG and MB-RWG basis functions.The new basis functions are grouped in two different blocks, including the MR basis functions defined at the lower levels and the gFs defined at the last level as where [T MR ] is the subblock describing the set of N MR basis functions defined at lower levels and T gF L is the subblock with the set of N gF L functions at the coarsest level, being N MR + N gF L = N , i.e., the total number of unknowns.In ( 14), [D] is a diagonal preconditioner whose elements are formally given by with T i the ith row of [T ].Note that although formally written like above, the diagonal preconditioner can be obtained directly from the diagonal of the near-field matrix in the context of a local MLFMA solution of the subdomain.
To further improve the conditioning of the new matrix system, an incomplete LU preconditioner is applied to the matrix block corresponding to the gFs [i.e., Z gF L in (14)].Considering all the above, the preconditioned problem can be written as with where [LU ] −1 is the inverse of the incomplete LU factorization of the matrix block corresponding to the gF L only.Finally, to speed up the evaluation of the global and local MVPs required in (10), a hybrid MPI/OpenMP parallel implementation of the MLFMA-FFT is applied in synergy with the DDM scheme [4].

IV. NUMERICAL RESULTS
A realistic numerical problem is presented to demonstrate the efficiency and versatility of the proposed approach to solve large problems with complex subsystems exhibiting local deep-multiscale features.All the simulations are performed with an AMD EPYC 7H12 64-CoreProcessor computing server with 4-TB memory.
Let us consider the characterization of the X-band system mounted on a morphed complex satellite, as shown in Fig. 1.The dimensions of the satellite are approximately 20 m length, 7-m beam, and 5 m height (800λ × 280λ × 200λ at the working wavelength, λ, at 12 GHz).The detailed model of this structure includes four different communication systems, as shown in Fig. 1.The X-band and Ku-band communication systems consist of four horn antennas supported by two external arms (shown in the top-left corner of Fig. 1) whose beam is directed by five reflectors.The Ka-band communication system consists of a complex array of horn antennas, which is shown in the bottom-right corner of Fig. 1.The model is completed with two solar panels and the main structure supporting the described elements.
The DDM is applied in conjunction with the proposed MR-LU preconditioner to get an accurate prediction of the radiation at 12 GHz Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Partition into subdomains of the satellite, highlighting some geometrical details.
of the X-band system 1 (see Fig. 1).A total of 60 267 853 unknowns (with 2160 MB-RWG) are required to mesh the entire geometry, which exhibits deeply multiscale features, including nonconformal discretization (shown in the top-right inset of Fig. 1), with a disparate mesh size ranging from λ/20 on smooth surfaces to λ/1350 close to the Ka-band antennas.This disparity in the mesh size results in poor conditioning, which in turn leads to lack of convergence unless proper simulation methodologies are applied.The excitation consists of a delta-gap voltage on the feed terminals of the antenna.
For the analysis of this challenging structure, the entire problem is partitioned into ten subdomains (also shown in Fig. 1).This partitioning is done according to its natural features: 1) two large subdomains corresponding to the complete solar panels, far enough from the excitation domain, where it is assumed a smooth variation of the electric currents and fast convergence; 2) four subdomains for the four aperture antennas belonging to the different communication systems on the satellite; 3) two subdomains for the support arms of the upper communication systems, including its horn antennas; 4) one subdomain for the main structure of the satellite; 5) one subdomain for the array antenna.Aiming for a global residual error below 10 −3 , the subdomains are solved using ten independent MLFMA-FFT solvers, with six levels for the support arms subdomains, five levels for the array antenna subdomain, and eight levels for the rest, setting a local residual error target of 10 −3 .The intrinsic complexity of this structure is exacerbated by the presence of the four communication systems, especially the large array with 90 horn antennas, shown in the right insets of Fig. 1.The strong mutual coupling between elements prevents, in this case, the use of 90 equal small MoM subdomains.Although this would allow to take advantage of the repetition pattern, relieving local calculations and memory consumption, it would be at the expense of the lack of convergence of the global iterative solver.Conversely, considering the entire array antenna as a single large subdomain is not without its drawbacks.The strong interaction between elements, its multiscale nature, and the use of nonconformal meshes will cause slow convergence, becoming a major bottleneck unless specific preconditioning techniques are applied.In this case, the proposed MR-LU technique is applied to speed up the solution of the subdomains.
We first analyze the residual error of the local iterative solution (at a given DDM iteration) of this horn array antenna subdomain, which is expected to give the worst convergence, as mentioned above.Fig. 2 shows the convergence performance (both in terms of iterations and in terms of wall-clock time) using the proposed  MR-LU preconditioner, compared to the diagonal, the SPAI [17], [18], and ILU [14] preconditioners built from the MLFMA nearfield matrix.It can be observed that the MR-LU preconditioner greatly reduces the solving time, allowing the effective solution of this challenging subdomain in the framework of the DDM iterative solution.
Next, we study the performance of solving the whole problem with DDM, by locally applying the three preconditioners mentioned in Fig. 2. Fig. 3 shows the comparison of the total time to solve by letting the local solvers to reach the target residual error of 10 −3 .A high reduction in solution time is observed with the application of the MR-LU preconditioner to local solutions, posing an accurate solution in less than 20 h.This is in contrast with the case of using diagonal local preconditioners, which would need more than 400 h, and the case of using ILU preconditioners, which would require more than 100 h to reach the final result.In the case of using the SPAI preconditioner, the convergence stagnates.
For a deeper comparison of the three configurations in terms of time and performance, we restrict the iterative solution of the local solvers by setting a residual error threshold of 5 × 10 −2 for the subdomains belonging to antennas.The results are also shown in Fig. 3, where, although a similar Krylov exterior convergence might be expected for all cases of local preconditioners, the poorer residual error exhibited by diagonal and ILU impairs the outer Generalized Minimal Residual Method (GMRES) convergence.This ruins the accuracy of the solution, turning the DDM with the proposed MR-LU preconditioner the only one of the considered approaches capable of adequately solving this challenging problem in a reasonable time.
For the sake of completeness, we include the convergence for the solution of the satellite applying the MR-LU preconditioner to the whole structure via the MLFMA solver (labeled as MLFMA-MR-LU); it exhibits a slower convergence.The computation times and peak memory for calculating the initialization for all the considered methods are shown in Table I, revealing that the DDM-MR-LU outperforms the other methods in both time and memory.
Finally, Fig. 4 shows the real part of the equivalent electric surface current distribution evaluated by applying DDM with the proposed MR-LU preconditioner.

V. CONCLUSION AND PERSPECTIVES
In this communication, we presented the efficient combination of the MR preconditioner with the DDM to obtain the accelerated solution of complex multiscale domains discretized with nonconformal meshes in the local stage.This strategy equips the local solvers with an efficient tool to address medium-sized complex subdomains exhibiting multiscale features, ensuring the optimal convergence of the DDM in real-life challenging problems.The efficiency and accuracy of the proposed methods were demonstrated through the solution of a realistic radiation numerical example in the field of space.
The next step of the work will be to extend the proposed method also to finite dielectric structures combined with metal parts.

Fig. 1 .
Fig. 1.Partition into subdomains of the satellite, highlighting some geometrical details.

Fig. 2 .
Fig. 2. Convergence performance of the inner GMRES when solving one outer GMRES iteration of the array of antennas domain of the satellite using the diagonal, the ILU, the SPAI, and the proposed MR-LU preconditioners in terms of (a) iterations and (b) wall-clock time.

Fig. 3 .
Fig. 3. Convergence performance of the outer GMRES when solving the EMC problem of the satellite using the diagonal, the ILU, the SPAI, and the proposed MR-LU in the local solver and using the MLFMA-MR-LU as global solver; ϵ is the threshold of the residual error in the local solvers.

Fig. 4 .
Fig. 4.Real part of the equivalent electric surface current distribution (dBµA/m) on satellite surfaces provided by the SIE-DDM approach.

TABLE I COMPUTATIONAL
TIMES AND PEAK MEMORY WHEN COMPUTING THE SETUP STEPS FOR THE DIFFERENT PRECONDITIONERS