{'fraction_non_alphanumeric': 0.04557145153512297, 'fraction_numerical': 0.02294647162835557, 'mean_word_length': 4.1596516690856316, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-performance Parallel Solver for Integral Equations of Electromagnetics Based on Galerkin Method

A new parallel solver for the volumetric integral equations (IE) of electrodynamics is presented. The solver is based on the Galerkin method which ensures the convergent numerical solution. The main features include: (i) the memory usage is 8 times lower, compared to analogous IE based algorithms, without additional restriction on the background media; (ii) accurate and stable method to compute matrix coefficients corresponding to the IE; (iii) high degree of parallelism. The solver's computational efficiency is shown on a problem of magnetotelluric sounding of the high conductivity contrast media. A good agreement with the results obtained with the second order finite element method is demonstrated. Due to effective approach to parallelization and distributed data storage the program exhibits perfect scalability on different hardware platforms.

geothermal or ore (high-conductive) reservoirs. Measured electrical and/or magnetic fields are further interpreted via the calculations for a given three-dimensional model of conductivity distribution. Maxwell's equation describing the EM field distribution can not be solved analytically in general case requiring numerical simulation. Large number of such simulations is required, and complex large scale models are invoked, Chave and Jones (2012).
The growing amount of data calls for the development of new numerical methods capable to deliver fast and accurate EM simulations and harness the computational power provided by modern high-performance multi-core and multi-node platforms.
The main difference in these approaches is in the discretization outcome. The model usually consists of a number of the non uniform three-dimensional anomalies embedded in the one-dimensional (layered) background media. While the FD and FE methods produce large sparse systems, Ernst and Gander (2012), the IE method results in the compact dense system matrices. The compactness is attained since the modeling region is confined only to the three-dimensional conductivity structures (anomalies) under the investigation, Raiche (1974); Weidelt (1975). Note that boundary conditions are satisfied by construction of the Green's functions. By contrast, in the FD and FE methods one has to discretize a volume much bigger in both lateral and vertical directions in order to enable the decay (or stabilization) of the EM field at the boundaries of the modeling domain, Grayver and Kolev (2015); Mulder (2006). Another distinction between the methods is in the condition number of matrices (it controls the stability of the solution). In FD and FE methods the condition number depends on the discretization and frequency, whereas in IE approach it does not, Pankratov et al (1995); Singer (1995).
The main focus in development and implementation of the numerical methods for the IE is the efficiency and performance for the large number of unknown parameters. The proposed new iterative numerical solver for the IE addresses this issue both mathematically (increasing the accuracy and stability of coefficient computations and reducing the memory usage) and computationally (allowing high degree of parallelism without memory loss at the nodes). The special class of integral equations is used: the integral equations with contracting kernel (CIE), Pankratov et al (1995); Singer (1995). The CIE were proved to have a unique solution and a wellconditioned (by construction) system matrix, Pankratov et al (1995); Singer (1995). The Galerkin method is used to solve CIE numerically, Delves and Mohamed (1985). In Kruglyakov (2011);Singer (2008) it was proved that the corresponding solution converges.
Two main issues in the IE numerical solution are the calculation of matrix coefficients with sufficient accuracy, Wannmaker (1991), and the storage of this dense matrix in some packed form, Avdeev et al (1997Avdeev et al ( , 2002. The Galerkin method with piece-wise constant basis allows to address both of these issues. Namely, the system matrix is decomposed into sums and products of diagonal and block-Toeplitz matrices. First, matrix coefficients (i.e., double volumetric integrals of the product of basis functions and CIE kernel) are analytically transformed into the one-dimensional convolution integrals. These convolution integrals are then computed by digital filtering approach. Keep in mind, that only weights of these filters are computed numerically, whereas the functions in the knots are computed analytically. The resulting system is solved using the Flexible GMRES, Saad (1993).
The proposed solver exhibits three main features. 1) The memory usage is 8 times lower, compared to the analogous IE solving algorithms. It is important to stress that the memory is saved for any background with an arbitrary number of layers. In contrast, in Kamm and Pedersen (2014) the memory reduction is achieved only in case of homogenous half-space as a background and for uniform vertical discretization. In Avdeev and Knizhnik (2009);Koyama et al (2008); Sun and Kuvshinov (2015) it is achieved at the expense of accuracy and performance. The idea behind the proposed solver is to combine the Galerkin method with the properties of EM field, namely Lorentz reciprocity. The matrix of the ensuing linear system can be separated in symmetric and antisymmetric submatricies. This reduces the memory requirements by 8 times, Sect. 3. 2) Efficient and accurate method for the computation of these matrices, Sect. 3. 3) The implementation with high degree of parallelism. The computational experiments performed with "Bluegene" and "Lomonosov" supercomputers from MSU, and high-performance computer (HPC) "Piz Daint" from Swiss National Supercomputing Center show that the solver makes the best usage of 128 to 2,048 nodes for calculation at a single frequency and a single source. The program exhibits perfect scalability. This paper is organized as follows. Section 2 is devoted to the overview of the CIE approach and the construction of the approximating system of linear equations. Section 3 addresses the reduction in memory requirements, the computation of matrix coefficients, and features of parallel implementation. In Sect. 4 the computational results for high (more than 3 · 10 4 ) conductivity contrast COMMEMI3D-3 model (Hursan and Zhdanov (2002);Varentsov et al (2000)) are compared with the corresponding results obtained using FE method by Grayver and Kolev (2015). The Appendices A to C provide the mathematical details of the presented method.

Overview
Assume that the EM fields are induced by the external electric currents J ext . Moreover, assume that the EM fields are time dependent as e −iωt , where ω is an angular frequency, i = √ −1 and magnetic permeability µ 0 is the same in whole space. Let σ (M), Re σ (M) ≥ 0 be a three-dimensional complex conductivity distribution in The solution of (1) is unique under the additional radiation conditions at infinity, Ward and Hohmann (1988).
Let Ω ⊂ R 3 be some bounded domain and σ (M) = σ b (z) for M(x, y, z) ∈ Ω and σ (M) = σ a (M) for M ∈ Ω (Fig. 1). Then for any M ∈ R 3 the fields E(M) and H(M) are expressed in terms of the integrals Here G E , G H are electrical and magnetic Green's tensors respectively, Dmitriev et al (2002);Pankratov et al (1995). The terms E N , H N are called the normal electric and magnetic fields, corresponding. They form the solution of the system with corresponding conditions at the infinity. Note, that G E , G H are indepent of the anomalous conductivity σ a . Let L 2 [Ω] be a Hilbert functional space of vector functions V with the following norm and dot product Suppose Re σ b (z) > 0 for M(x, y, z) ∈ Ω (the typical EM sounding situation), then the operator G m E is defined as where G E is an integral operator from the first equation in (2). The operator G m E is a contracting operator in L 2 [Ω], Pankratov et al (1995);Singer (1995). Using (2) and (5) one obtains the CIE for E where I is the identity operator and σ b means complex conjugation of σ b .
Let W N be a linear span of the vector functions W n , W n = (W n x ,W n y ,W n z ), n x , n y , n z = 1 . . . N and P N be a projection operator from where γ = x, y, z. Note that P N = 1. Applying P N to the first equation in (6) one obtains the operator equation in W N Since b a < 1, G m E is a contracting operator and P N = 1, it can be easily shown that (9) has a unique solution W in W N , Kruglyakov (2011);Singer (2008 Using the definition of W N the components of W = (W x ,W y ,W z ) can be expressed as Using (8) to (10) and taking into account that σ a , σ b are piecewise functions, one obtains the following system of linear equations for the coefficients U n = (U x n ,U y n ,U z n ), where Note that K m n , B m n ,Î,γ m are 3 × 3 matrices,Î is an identity matrix,γ m is a diagonal matrix. The system (11) has a unique solution, Kruglyakov (2011);Singer (2008).
Using the solution U n of system (11) one can approximate E(M) and H(M) for any point M ∈ R 3 with Relations (13) are first order approximations of d in C for M ∈ Ω, Kruglyakov (2011), and result in fast and relatively simple computations. The main challenge is to calculate matrix coefficients and solve the system (11).

Memory Requirements
The main challenge of the integral equation approach is in solving of the system of linear equations with dense matrices (11). The storage of these matrices in RAM is also problematic. The standard approach, Avdeev et al (1997), is to use the property For the implementation purposes consider now Ω ⊂ R 3 to be a rectangular domain. As before Ω is divided in N = N x N y N z rectangular subdomains Ω n , n = 1, . . . , N, where N x , N y , N z are the number of subdomains in X,Y, Z directions respectively. Suppose also that each Ω n has the same size h x × h y in XY plane. Then

is also a block Toeplitz matrix and is induced by the block vector
Here Q αβ are the matrices of the order N z , α, β = x, y, Let A be a matrix corresponding to the system of linear equations (11). Then where S, R 1 , R 2 are the diagonal matrices; B = B m n is the block Toeplitz matrix described above.
In view of (17) it follows that only 36 · N x N y N 2 z · 16 + O(N x N y N z ) bytes are required to store matrix A in double precision. Using the equivalence G E xy = G E yx this requirement can be reduced to 32 · N x N y N 2 z · 16 + O(N x N y N z ) bytes as in Avdeev et al (1997). This memory requirement can be reduced in 8 times by virtue of the following Lemmas.
is an electrical Green's tensor of any layered media, then it possesses symmetric and antisymmetric properties in Cartesian coordinates along the vertical dimension is an electrical Green's tensor of any layered media, then it possesses symmetric and antisymmetric properties in Cartesian coordinates along the lateral dimensions These Lemmas are trivial corollaries from Lorentz reciprocity, Ward and Hohmann (1988) and formulas for the Green's tensor components (See Appendix A). Relations (15) and (18) give where T indicates a matrix transpose. Therefore, one needs to store only Q xz , Q yz and upper diagonal parts of Q xx , Q xy , Q yy , Q zz . Moreover the values Q(i, j) can be stored only for i = 0, . . . , N y − 1, j = 0, . . . , N x − 1, since (19) allows to obtain these values for negative i or j from suitable symmetric/antisymmetric properties.
Thus only 2·N x N y N z ·(2N z +1)·16 bytes are required to store B m n which is 8 times less than the memory requirements in Avdeev et al (1997Avdeev et al ( , 2002; Hursan and Zhdanov (2002). It is worth to stress again, that this is valid for any background layered media and without the conditions on the subdomains to be of the same vertical sizes.

Matrix Coefficients Computation
The next computational challenge of the Galerkin approach is the evaluation of the coefficients B m n , n, m = 1 . . . N, that is the double volumetric integrals of the G E in the RHS of (15), with desired accuracy. The components of G E are the improper integrals containing the Bessel functions, Appendix A. The integration in vertical direction is performed analytically using the fundamental function of layered media approach from Dmitriev et al (2002), Appendix B. The main problem, however, is the integration over the horizontal domains. In this case one needs to compute the fifth-order integrals over the fast-oscillating functions.
The integrals in (15) are double volumetric ones, thus they have only weak singularity. Therefore, one can change the order of integration and make an appropriate substitution and convert the fifth-order integral to a convolution with the specific kernel. Following the standard approach of convolution calculation the spectrum of this kernel is computed and the digital filter is constructed, Appendix C.
It is important to emphasize that both the knots and the weights in the obtained filter significantly depend on the integration domains. On the contrary, the integration over different horizontal domains is completely data independent. This is used in parallel algorithm. The computational experiments demonstrate (Sect. 4) that the used filters provide suitable accuracy even for the models with high conductivity contrast.

Parallel Implementation
The most essential part of any iterative method for solving a system of linear equations is the matrix-vector multiplication. Since matrix B is a block Toeplitz matrix, one can use the two-dimensional Fast Fourier Transform (FFT) to speed up this operation, Avdeev et al (1997). Therefore, instead of matrices Q(i, j) the discrete Fourier transformations Q(i, j) are stored. This requires the same amount of memory since the discrete Fourier transform preserves the symmetric/antisymmetric properties of data.
The multiplication of block Toeplitz matrix B on some vector V ∈ W N is performed via the following three-step algorithm: 1. Compute 3N z forward FFT of vector V; 2. Compute 36N x N y algebraic matrix-vector multiplications of order N z to obtain vectorṼ; 3. Compute 3N z backward FFT of vectorṼ.

The multiplications in
Step 2 are further divided into 4N x N y groups which are mutually data independent. This allows to implement the special scheme of distributed data storage and a parallel algorithm of IE solver, described below.
For simplicity, consider 2N y nodes and assume that N x is even. The distributed storage of matrix is organized in a special way: the half of block-vector Q(n, j), j = 0 . . . N x /2 − 1 is stored at nth node, n = 0 . . . N y − 1, while Q(n, j), j = N x /2 . . . N x is stored at node n + N y , n = 0 . . . N y − 1, Table 1. (i) The coefficients of matrices B, S, R 1 , R 2 stored at different nodes are computed simultaneously and completely data independent; (ii) The iterative method is executed using the authors' distributed implementation of FGMRES Saad (1993), inspired by Fraysse et al (2003); (iii) The distributed two-dimensional Fourier transform is computed via the authors' implementation using FFTW3 library Frigo and Johnson (2005) for local FFT; (iv) The calculation of the local algebraic matrix-vector multiplication is processed by using OpenBLAS library Wang et al (2013) ; (v) For all the stages of the computational process the hybrid MPI+OpenMP scheme is used.
To demonstrate the scalability of the implemented parallelization the COMMEMI3D-3 model is used with N x = 176, N y = 224, N z = 118, that is with cubic subdomains with 25 m edges, Sect. 4. The computational experiments performed at "Bluegene/P", HPC "Lomonosov" (MSU) and Piz Daint (Swiss National Supercomputing Center) showed good speed increment depending on the number of processes (Fig. 2). Matrix calculation time includes time of FFT calculation of B m n . The solid black line means ideal linear speed up. Note, that for such high-contrast model matrix calculation time (crosses) is small enough compared to solving of the system of linear equations (circles). One can see that the scalability is close to a linear.

High Conductivity Contrast Modeling
The accurate computation of the EM field in a high conductivity contrast media is one of the most complex problems of EM modeling due to strong codependency between conductivity contrast and matrix condition number, Pankratov et al (1995);Pankratov and Kuvshinov (2016);Singer (1995). The conductivity contrast means the ratio between the real parts of anomalous conductivity Re σ a (M) and background conductivity Re σ b (M) at the same point M.
The high conductivity contrast COMMEMI3D-3 model, Hursan and Zhdanov (2002); Varentsov et al (2000), is used as one of the test models for the presented solver. This model schematically describes the conductivity distribution typical for  Table 2: The coordinates of the opposite corners (x 1 , y 1 , z 1 ), (x 2 , y 2 , z 2 ) in km and conductivities σ (S/m) of COMMEMI3D-3 blocks the ore exploration by the audio-magnetotelluric sounding. Following magnetotelluric (i.e., low-frequency) sounding tradition in the rest of this section the conductivity is a real-valued function. The COMMEMI3D-3 model consists of seven rectangular blocks placed in a layered media and oriented along coordinate axes. Their conductivities σ (in S/m) and positions (coordinates of the opposite corners in km) are listed in Table 2 and depicted in Fig. 3. The layered background of the model consists of the upper halfspace z < 0 (air) with conductivity of 0 S/m, two layers with the conductivity of 10 −3 , 10 −4 S/m, and a lower halfspace with conductvity of 0.1 S/m. The thickness of the first and the second layers is 1 and 6.5 km respectively. One can see that maximum conductivity contrast is 10 4 in the first layer and 3.3 · 10 4 in the second one.
The modeling of magnetotelluric sounding was performed for various discretizations (cubic subdomains of different sizes) and periods T = 2π/ω, and was compared with the results from modern FE solver by Grayver and Kolev (2015). Figures 4 and 5 represent, correspondingly, the apparent resistivities ρ xy at profile x = 1.9 km and ρ yx at profile y = 3.83 km for the period 1 s. One can see that the agreement with FE (magenta circles) is good even for rather coarse anomaly discretization (black curve). The exception is the area [3.5, 4.5] km on the profile y = 3.83 km (Fig. 5) above the high-conductivity block. This is amended by taking finer discretization (azure curve).  Figure 6 shows the apparent conductivity ρ yx at site (3.975, 3.83) km, that is above the high-conductivity block, at different periods. One can see that the finer discretization is needed only for the periods 10 −1 , 10 1 s. The reason behind this effect is the drastic change in the electric field inside of the compact high-conductivity block that does not allow to use the piecewise approximation on coarse discretization. The solver by Grayver and Kolev (2015) uses the second order polynomials which are very effective in such situations. At the same time, the coarse discretization can be efficiently used for smaller and larger periods. It is worth emphasizing that this concerns only the area above the high-conductivity block, while at the point (1.9, 1.7) km, Fig. 7 shows good correspondence for all periods.
Figures 8 to 11 demonstrate area distribution of the apparent conductivity and impedance phases for 1 s period. One can see that the variation in apparent conductivity is of the four order of magnitude with very drastic transition. The phase change of the impedance is 10 degrees, and the transition is again drastic. It is worth reminding, that one of the peculiarities of IE method is quite weak dependence of computational costs on the number of sites where the field is computed. This allowed to obtain the maps with such drastic transitions without the additional computational costs.

Conclusion
The presented solver named "Gnu Integral Equation Modeling in ElectroMagnetic Geophysics" (GIEM2G) shows impressive performance in terms of both memory requirements and accuracy. The memory requirements are 8 times lower compared to other volumetric IE solvers Avdeev et al (1997), Hursan and Zhdanov (2002). It is achieved for any layered background and non uniform discretization in vertical direction. In this way the average-scale modeling (up to 3 · 10 6 subdomains) can be efficiently done using laptops. The parallelization scheme allows to use HPC with hundreds and thousands of nodes for large-scale modeling (up to 10 9 subdomains).
The computational efficiency of the method is demonstrated on high-conductive contrast (3.3 · 10 4 ) model COMMEMI3D-3. To the best of the authors knowledge, it is the first time that such high-contrast complex model provides comparable results for such different methods as FE and IE. In addition to the efficient usage of HPC the proposed IE method relies on the new technique to calculate the matrix coefficients. It is based on the analytical integration in vertical direction and completely new scheme to compute the integrals in horizontal direction. It is worth mentioning that the proposed scheme of analytical integration is robust in terms of machine precision and needs only O(N z ) computations of complex exponents.   From (21) and (22) one gets Let , where S n , S m are horizontal rectangular domains, and let σ That is the subdomains do not intersect the boundaries of the layers. Taking into account (12), (21) and (24) one can see that B m n is expressed in terms of double volumetric integrals with weak integrable singularity, so the order of integration can be changed. Then using (23) and (24) where The function U γ (z, z * , λ ) is defined as a unique solution of the problem The following procedure is performed to obtain an explicit expression for U γ that allows analytical integration. Let l 0 = 0, l N lay = 0, l n = d n+1 − d n , n = 1 . . . N lay − 1. Define p γ m , q γ m , m = 0 . . . N lay by the recurrent expressions Let d 0 = d 1 , d N lay+1 = d N lay and let points z r , z s belong to r and s layers respectively, 0 ≤ r, s ≤ N lay . Then using (29) one gets r,s p γ r e 2η r d r e −η r z r + e η r z r e −η s z s + q γ s e −2η s d s+1 e η s z s for z r ≤ z s ; A γ s,r p γ s e 2η s d s e −η s z s + e η s z s e −η r z r + q γ r e −2η r d r+1 e η r z r for z r > z s , where A γ r,s = Q γ r × Q γ r+1 × · · · × Q γ s−1 A γ s,s , for r < s, γ n e −2η n l n , for r = s = n, n = 0 . . . N lay , To check (30) one can explicitly substitute (30) in (28) taking into account (29) and (31).
In view of (30) one can see that integrals in (27) (i.e., the integrals over U γ (z, z * , λ ) and its partial derivatives) can be integrated analytically with respect to z, z * over any domains that do not intersect the layer boundaries. However, the rounding errors arising in addition and multiplication of very small or large quantities make the formula (30) impractical for λ 1. Instead the following formula is used r,s p γ r e −(η r z r +η s z s −2η r d r ) + q γ s e −(2η s d s+1 −(η r z r +η s z s )) + e −(η s z s −η r z r ) + p γ r q γ s e −(2(η s d s+1 −η r d r )−(η s z s −η r z r )) for z r ≤ z s ; A γ s,r p γ s e −(η s z s +η r z r −2η s d s ) + q γ r e −(2η r d r+1 −(η s z s +η r z r )) + e −(η r z r −η r z r ) + p γ s q γ r e −(2(η r d r+1 −η s d s )−(η r z r −η s z s )) for z r > z s .
Formula (32) overcomes the aforementioned problem, since the real parts of all the exponents powers are negative. The consequent calculations provide accurate and robust results for any 0 < λ < ∞.
Consider N z subdomains in the discretization in vertical direction. To obtain the matrix B m n for the system (15) one needs to compute O N 2 z complex exponents in (32). An algorithm requiring only O (N z ) complex exponents calculations is developed to speed up the integration procedure.
For fixed p, q, ϕ the integral in (37) is the convolution integral with kernel Φ. Note that for different values of α and β the kernels can be obtained similarly.
The main advantage of using the convolution integrals is that their computation does not require the explicit calculation of kernel Φ. Consider the input function v(t) and the output function u(s) such that Then for any g(t) one gets The tradeoff between the accuracy of (40) and the computational time is achieved by the particular selection of M, N 1 , N 2 , l, ξ and functions u and v.
From (36) and (40) the approximation formulas for (35) can be obtained For the given input function v(t) = 8e −t 2 t 5 − 4t 3 + 2t the output functions for different kernels can be expressed analytically by Gaussian and error functions. Inspired by Anderson (1979) the parameters used are l = 0.2, ξ = 0.0964, M = 512, N 1 = −250, N 2 = 200. In computational experiments these parameters provided appropriate accuracy in calculation of B m n .