Iterative solution methods for 3D controlled-source electromagnetic forward modelling of geophysical exploration scenarios

We develop an efficient and robust iterative framework suitable for solving the linear system of equations resulting from the spectral element discretisation of the curl-curl equation of the total electric field encountered in geophysical controlled-source electromagnetic applications. We use the real-valued equivalent form of the original complex-valued system and solve this arising real-valued two-by-two block system (outer system) using the generalised conjugate residual method preconditioned with a highly efficient block-based PREconditioner for Square Blocks (PRESB). Applying this preconditioner equates to solving two smaller inner symmetric systems which are either solved using a direct solver or iterative methods, namely the generalised conjugate residual or the flexible generalised minimal residual methods preconditioned with the multigrid-based auxiliary-space preconditioner AMS. Our numerical experiments demonstrate the robustness of the outer solver with respect to spatially variable material parameters, for a wide frequency range of five orders of magnitude (0.1-10’000 Hz), with respect to the number of degrees of freedom, and for stretched structured and unstructured as well as locally refined meshes. For all the models considered, the outer solver reaches convergence in a small (typically < 20) number of iterations. Further, our numerical tests clearly show that solving the two inner systems iteratively using the indicated preconditioned iterative methods is computationally beneficial in terms of memory requirement and time spent as compared to a direct solver. On top of that, our iterative framework works for large-scale problems where direct solvers applied to the original complex-valued systems succumb due to their excessive memory consumption, thus making the iterative framework better suited for large-scale 3D problems. Comparison to a similar iterative framework based on a block-diagonal and the auxiliary-space preconditioners reveals that the PRESB preconditioner requires slightly fewer iterations to converge yielding a certain gain in time spent to obtain the solution of the two-by-two block system.


Introduction
Geophysical electromagnetic (EM) methods are crucial tools in mapping the Earth's subsurface in terms of its electric conductivity distribution. These methods are particularly employed for geophysical surveys targeting hydrocarbon and mineral resources (cf. e.g. [51,86]), or for investigations concerned with geothermal, hydrogeological or environmental applications (cf. e.g. [65,82,87]). Such surveys have recorded large amounts of three-dimensional (3D) EM data in vastly different geological settings. The interpretation of such acquired data sets involves solving an inverse problem to obtain a subsurface model that best fits the recorded data. Although the inverse problem entails considerable computational burden itself, the repeated solution of the 3D EM forward problem within the inversion framework constitutes one of the core factors for the high computational cost of the inversion. This is further aggravated considering that realistic surveys compromise large 3D volumes with many receivers and possibly multiple sources demanding large 3D meshes to accurately represent realistic subsurface structures and topographic variations yielding forward problems of large proportions, often including tens to hundreds of millions degrees of freedom. In light of this, solving the forward problem efficiently and effectively is crucial to limit the computational effort, in particular the time and memory requirements, for the forward problem itself as well as for inversions of realistic industrial 3D data sets.
The 3D EM forward simulation consists of solving the governing partial differential equation (PDE) derived from of Maxwell's equations using discretisation. In EM geophysics, a number of different approaches have conventionally been employed in the past few decades to discretise the relevant continuous PDEs. Such methods include the finite difference (FD), finite element (FE) and finite volume (FV) methods. Regardless of the used approach, the corresponding discrete equations are expressed in terms of a linear algebraic system of equations of the form Ax = b, where A is a large-scale sparse nonsingular matrix. There exists numerous numerical methods for solving the resulting linear system, which are of two general types -direct and iterative, both with their inherent advantages and drawbacks.
It is well known that obtaining accurate and efficient solutions of the algebraic system of equations arising from Maxwell's equations is a computationally difficult task. This is due to several factors which impact the condition number of the system matrix A negatively and can make it very ill-conditioned : (1) the influence of the large kernel of the curl operator ∇× present in Maxwell's equations [80,81]; (2) large electric conductivity contrasts, in particular across the air-Earth interface where the contrast amounts to four orders of magnitude or more; (3) the requirement of non-uniform and locally refined meshes; (4) using high-order polynomials p for the space approximation not only yields a more densely populated system matrix but also a more ill conditioned one. In fact, the condition number of the system matrix in 3D increases with O(p 8 ) for the finite element method as shown by [1].
Direct methods are based on an exact lower-upper triangular (LU) factorisation of the matrix A or its Cholesky decomposition LL T when A is symmetric and positive definite. Examples of efficient sparse direct solvers that are widely employed in the geophysical EM community include a distributed-memory multifrontal massively parallel sparse direct solver called MUMPS [3,4] and a shared-memory parallel direct solver known as PARDISO [77]. Direct methods are robust and accurate, even for systems with illconditioned matrices. In addition, for the solution of linear systems with multiple right-hand sides as in multi-source surveys in controlled source electromagnetics (CSEM), direct solvers only require a single factorisation of the system matrix (see e.g. [32,66,85]). On the other hand, the memory demands and time requirements grow as O(N 2 ) and O(N 4/3 ) in 3D, respectively, where N is the number of unknowns and can limit the size of feasible problems depending on the available computer platform. In regard to (parallel) scalability, direct solvers are considered nonoptimal on distributed-memory platforms due to a low computation-to-synchronisation ratio [43,61]. In practice, this means that, due to the additional overhead caused by communication and synchronisation between the nodes, beyond a certain number of parallel processes direct computations cannot be sped up by adding more computing nodes/cores (see e.g., [32,68]). Nevertheless, owing to successive advances in direct solution algorithms [5,78], as well as the increasing availability and power of highperformance computing facilities [60], direct solvers are often used when solving the EM forward problem see (e.g. [29,32,35,68,71,72,81,85], for details).
Iterative methods, e.g. Krylov subspace iteration methods, include a broad range of solution methods that use successive approximations to acquire progressively more accurate solutions of the considered linear system of equations at each step until reaching convergence. The main advantage of iterative methods is that they require less computational resources than direct methods in terms of time and memory, and are thus better suited to handle large-scale problems. The major computational building blocks of iterative solvers include cheap matrix-vector products and vector operations which can be easily parallelised making them more scalable. However, for ill-conditioned systems, iterative methods, if straightforwardly applied, can exhibit slow convergence or even diverge. Therefore iterative solvers require the usage of an efficient preconditioner. Preconditioners transform the original system of equations into an equivalent one that has a (preconditioned) matrix with a more favourable eigenvalue spectrum. This ought to improve the convergence rate of iterative methods whilst sufficiently offsetting the additional computational burden of building and applying the preconditioner. In this way, preconditioned iterative solvers become viable and cost-effective solution methods, in particular from the perspective of large-scale 3D problems. The task to construct a numerically effective and computationally efficient preconditioner is far from trivial and in many cases the preconditioner might be problemspecific.
Examples of iterative methods and corresponding preconditioners include the Biconjugate gradient (BICG) or the Quasi-minimal residual (QMR) method with different preconditioners such as diagonal or Jacobi scaling [2,28], the incomplete Cholesky decomposition [83], or the incomplete LU factorisation [89]. Other commonly used iterative solvers are the Generalised minmial residual (GMRES) and the Biconjugate gradient stabilised (BICGStab) methods which, for instance, have been used by [39] and [49] preconditioned by incomplete LU decomposition. [67] investigate the convergence behaviour of the aforementioned iterative methods preconditioned with a sparse approximation of the inverse of the system matrix based on the minimisation of a suitable Frobenius norm. In addition and in order to solve problems with multiple right-hand sides (RHS), [66] present block GMRES and QMR methods [76] using either an incomplete LU, an approximate inverse or a simplified geometric multigrid preconditioner.
The class of geometric Multigrid (MG) and algebraic Multigrid (AMG) methods has gained popularity, because in many cases they are optimal regarding numerical efficiency, they converge in a number of iterations, that is independent of the number of degrees of freedom, and they also have optimal (or nearly optimal) computational complexity, namely, linear with respect to the number of degrees of freedom. MG methods make use of changing discretisation on a sequence of successively coarser grids. These techniques are based on the observation that although most relaxation methods, e.g. Jacobi and Gauss-Seidel algorithms, may in general converge slowly, they attenuate effectively and very rapidly the oscillatory (high-frequency) modes of the errors of any given mesh. However, standard relaxation methods hardly damp the smooth (low-frequency) components of the error within the few first iterations. Thus, after removing the oscillatory error components, the convergence of these basic iterative methods slows down as the smooth modes are eliminated slowly. However, the smooth modes on a given discretisation mesh naturally become oscillatory modes on a coarse grid, which is the origin of the idea of moving to a coarser mesh discretisation in order to remove the corresponding error mode once the relaxation methods start stalling. The process of moving to coarser grids can be repeated recursively, in this way attenuating effectively and rapidly the oscillatory modes on the different levels, thus reducing the global error quickly. In contrast to MG methods, AMG methods achieve the same effect however, without using any discretisation meshes but creating an analogous algebraic structure based on an aggregation technique. Therefore, AMG methods are more generally applicable than MG methods and are their preferred alternative. Multigrid methods can be used as stand-alone solution methods but more often are used as preconditioners. These methods are considered scalable [31,94] and can have linear complexity O(N) in computational and memory load for various boundary value problems [36,88] making them a viable option for large-scale simulations. Despite this and in comparison to direct and iterative methods, multigrid methods have not been widely used in geoelectromagnetic applications, which might be due to their involved numerical implementation and the fact that generic multigrid schemes fail for Maxwell's equations due to the large kernel of the curl operator [26,80,81]. Since the early 2000s, renewed development in multigrid methods, in particular for the curl-curl operator [46][47][48]70], as well as the availability of several software libraries such as hypre [38] and ML [40], have improved the applicability of MG methods to geophysical EM problems. For example, [6] apply a geometric multigrid preconditioner based on Dendy's black box multigrid solver [33]. Mulder [58] has developed a geometric multigrid method which is either used directly or as a preconditioner and has reported good convergence rates for uniform grid spacing, but inadequate rates for substantially stretched grids. More recently, [52] have implemented an AMG scheme as a preconditioner for Krylov subspace iteration solvers that is shown to provide mesh-independent rate of convergence.
Often MG and AMG methods are used in various preconditioning techniques as a block-solver. An example of this for application with conductors can be found in [30] in the following context. A complex-valued system is rewritten as a real two-by-two block system, preconditioned by a block-diagonal preconditioner. The two linear systems of the diagonal blocks are solved by an iterative method with the auxiliary-space technique [48,54]. The main idea consists of approximately inverting the curl-curl operator by casting the problem into subspaces in which it can be solved efficiently by AMG methods. This preconditioning framework introduced in [30] has also been adopted to EM modelling problems (cf e.g. [25,44,45,69]).
Two-by-two block systems are encountered in numerical applications in various fields, e.g. discrete Navier-Stokes problems, eddy current problems, linear elasticity, optimisation problems constrained by PDEs etc. As mentioned, these also arise when transforming a given complex-valued system to its real-valued equivalent. As the occurrence of block systems is quite common much research has been carried out and hence efficient preconditioners for block-structured systems can be found in the literature. Here, we restrict our selection to a few robust and efficient techniques for block-structured systems. These include block-diagonal and block-triangular preconditioners [12,56,84] and the preconditioned modified Hermitian -skew Hermitian splitting (PMHSS) iteration method [17][18][19]. Another preconditioning technique presented in various publications (cf. e.g., [9,10,15,27]) is referred to as the PREconditioning for Square Blocks (PRESB).
The objective in this work is to present a robust and efficient iterative solution framework for the algebraic system of equations arising from discretising problems in 3D frequency-domain controlled-source electromagnetic modelling. For the solution of the system, we consider the PRESB preconditioning technique [10] to be used in a Krylov subspace iteration method. An action of the PRESB preconditioner amounts to solving two additional linear systems, which we solve by an inner iterative method preconditioned by the auxiliary-space preconditioner [48,54]. We emphasise that PRESB, although resembling the block-diagonal preconditioner presented in some other related work (e.g. [25,44,45,69]), is different. In particular, applying these preconditioner results in slightly different computational procedures although the core of both schemes consists in solving two linear systems. Further and at a later point in this work, we show that the PRESB and block-diagonal preconditioners yield different eigenvalue bounds for the preconditioned block systems. To the best of our knowledge, so far this preconditioner has not been employed in geo-EM related applications. We also stress that we use a spectral element discretisation method, whereas previous works use finite element methods. In contrast to [44,45] and [69], we do not disregard displacement currents which complicates the system due to the potential negative shift in matrix elements resulting from the displacement current term in some parts of the modelling domain [25].
The remainder of this paper is structured as follows. Section 2 describes the mathematical model, its discretisation and the properties of the arising algebraic system of equations. In Section 2.1, we present the formulation of the continuous problem. We then derive the dimensionless form of the problem at hand in Section 2.1.1 and briefly describe its discretisation using the spectral element method in Section 2.2. Section 3 describes the solution procedure we choose and some possible block preconditioning techniques for two-by-two block system, such as the blockdiagonal, the Schur complement-based and the PRESB preconditioners in Sections 3.2, 3.3 and 3.4, respectively. When introducing the PRESB preconditioner, we also present the general block iterative framework. Section 3.4.1 briefly describes the auxiliary-space Maxwell preconditioner. Numerical results for two test problems illustrating the robustness of the solver are presented in Section 4. We compare the performance of the PRESB preconditioner to the diagonal preconditioner employed in [44,45] and [69]. Some concluding remarks are given in Section 5.

Formulation of the continuous problem
Starting from Maxwell's equations in frequency domain the governing partial differential equation for the total electric field in CSEM applications is derived by eliminating the magnetic field and reads as follows, where is the bounded model domain, E denotes the total electric field, ω is angular frequency, μ is magnetic permeability, σ is electric conductivity, ε is dielectric permittivity and J s is the known imposed source current density. In electromagnetic geophysics, all three material properties μ, σ and ε are assumed to be spatially varying as well as piecewise smooth. However, we note that both the magnetic permeability and dielectric permittivity occupy a relative narrow range. The parameter μ can be safely assumed to lie within 1μ 0 to 3μ 0 , where μ 0 is the magnetic permeability of vacuum μ 0 = 4π ·10 −7 H/m. The parameter ε ranges from 1ε 0 to 81ε 0 (relative dielectric permittivity of water) with ε 0 being the dielectric permittivity of free space ε 0 = 8.854 · 10 −12 F/m. The remaining physical property σ , on the other hand, is strongly material-dependent and can range from 10 2 S/m for highly conductive massive sulfide ore deposits to 10 −6 S/m for very resistive metamorphic and igneous rocks. In addition, most applications in EM geophysics set the air conductivity on the order of 10 −7 to 10 −10 S/m to ensure that it is sufficiently small but larger than zero.
The range of frequencies employed in CSEM depends on the objective of the survey as well as the system available and its specifications. However, the frequency range influences the size of the computational domain. In particular, to limit computational boundary effects, the domain has to be chosen adequately. Typically, the planewave skin depth criterion which is a function of the background conductivity σ b , permeability μ b , permittivity ε b and the employed source frequency f can be used as a guideline to set appropriate domain boundaries. In order to adequately attenuate the electric field at the domain boundaries, the computational domain is extended up to a distance of approximately five plane-wave skin depths from the area of interest. It is further noted that the frequencyand material-dependent wavelength influences the choice of the element size as the grid spacing needs to be sufficiently fine compared to the wavelength in order to ensure accurate numerical computations [57]. Grid spacing is typically chosen as a fraction of the skin depth and based on element order. For first-order approximation, guidelines on grid spacing relevant for finite difference schemes and for finite element modelling in geophysical electromagnetics have been presented in [92] and [91] respectively.
In order to guarantee the uniqueness of the electric field E, Eq. 2a is subjected to homogeneous Dirichlet boundary conditions on the domain boundary ∂ , namely, (2b)

Dimensionless formulation of the curl-curl equation
Given the vastly varying orders of magnitudes of some of the model parameters and the large scale of the space domain, it is preferable to introduce dimensionless quantities for the governing Eq. 2a. The dimensionless forms of the equations serve to properly scale the problem and avoid ill-conditioning due to scaling -the problem becomes independent of measurement units and facilitates better understanding of the nature of the equations without involving physical parameters (see e.g. [22,24], for more details).
To obtain a dimensionless form we introduce the following new dimensionless variables, denoted by . and the corresponding suitable scaling factors, denoted by subscript f , namely, where x is a three-dimensional space vector. The curl operator ∇×, being a spatial derivative, transforms in a similar way, We now replace the dimensional quantities in Eq. 2a with their respective dimensionless equivalents yielding the dimensionless curl-curl equation and collecting the scaling factors for each term in dimensionless coefficients results in where it can readily be seen that the coefficients α, β and γ are indeed dimensionless and defined as follows The individual scaling factors have to be chosen in a suitable manner. Spatial dimensions are rescaled by L f which is taken to be the longest spatial dimension of the model. The latter ensures that the characteristic size of the space discretisation parameter, usually denoted by h, is less than 1. Further, the model conductivities are regularised by choosing σ f to be equal to the largest conductivity of the conductivity distribution. Similarly, the dielectric permittivities ε and the magnetic permeabilities μ are rescaled with their respective largest value. This means that the largest dimensionless conductivity σ, permittivity ε and permeability μ values are equal to one. The remaining scaling coefficients ω f , E f and J s,f are set to one. At this point we also note that, after solving the dimensionless system, to obtain the unscaled electric field we have to rescale the solution vector by 1 To simplify the notations, in the sequel we drop and, unless stated otherwise, we work with the dimensionless parameters and the scaled variables.

Discretisation approach
In this work, we use a particular discretisation approach, referred to as the Spectral Element Discretisation, described in detail in [93]. The main reason for using the spectral element method instead of the more traditional finite element method is to combine and harness advantages of the spectral method in form of high accuracy and favourable convergence properties [42,62,73] while retaining the geometrical flexibility of the finite element method (see e.g. [62,73,93,95]). For completeness, we include a short description below.
Following Galerkin's method, by taking the L 2 -inner product f, g L 2 := f · g of the dimensionless curlcurl equation for the electric field Eq. 4 with a vector test function , integrating over and using integration by parts, one obtains the weak form of Eq. 4 (cf. [57]): for all ∈ H 0 (curl, ), and where H 0 (curl, ) is defined as The space H(curl, ) denotes the space of functions defined on the bounded domain , which are in the square integrable space (L 2 ( )) 3 and such that the curl is well defined in (L 2 ( )) 3 The spectral element discretsiation of the weak form Eq. 5 of the boundary value problem Eq. 2a starts by subdividing the domain into a set of non-overlapping hexahedra. In addition, one needs to replace the function space H 0 (curl, ) by a finite dimensional approximation space consisting of curl-conforming vector functions which are element-wise polynomial functions of arbitrary degree within each element and non-zero for all but a few hexahedra [81]. We note, that due to the specific functional space, the tangential components of the vector basis functions are continuous across the interface of two neighbouring elements, whereas their normal components may be discontinuous across the same element boundary [59]. Letting H 0,h be such an approximate subspace to H 0 (curl, ) defined on the mesh, one can construct a finite set of basis functions that spans this subspace. This then yields the discrete formulation to find the approximate solution E h ∈ H 0,h such that for all ∈ H 0,h , and where the more compact notation of the inner product is used. Expanding the solution as to be fulfilled by any l, or in matrix-vector notation where C A is the sparse symmetric complex-valued system matrix, E h is the solution vector collecting the basis function expansion coefficients E k which are also referred to as degrees of freedom (DOF) and n dof is the total number of these degrees of freedom. The right-hand side b contains the source term entries. Further, K denotes the sparse symmetric positive semi-definite stiffness matrix and M σ and M ε are the sparse symmetric positive definite mass matrices defined by We note that the basis functions are real-valued and, hence, so are the matrices. As to the choice of the basis functions for the spectral element discretisation, we refer to [93] for details.

Iterative solution and preconditioning approaches
As already noted, the matrix (10) is complex symmetric and of very large dimension. Due to the problem parameters, it is also very ill conditioned. Therefore, we need to solve the system with some iterative method and use an efficient preconditioner. We aim at achieving a preconditioned iterative solver, which is robust with respect to the problem parameters ε, σ, μ, ω and also with respect to the discretisation parameter h.
The system (10) can be solved in two ways. As a first alternative, it can be solved using complex arithmetics, utilising an appropriate Krylov subspace iterative method, such as GMRES [76], BICG [90] and some others. However, the question how to choose a good preconditioner remains open. Standard choices, such as incomplete LU techniques, are known not to be robust. Therefore, we choose the second alternative, i.e. to rewrite the system in real form of twice larger size.
As is well known, any n × n complex system of the form (A + iB)(x + iy) = (u + iw) can be transformed into a 2n × 2n real system in several equivalent ways. A straightforward computation shows some of the equivalent forms of the arising two-by-two block real systems, Which form to use may depend on the properties of the matrices A and B. In our case the block M σ is symmetric and positive definite and we recast Eq. 10 into a realequivalent form by splitting the electric field into real and imaginary parts E h = E Rh + iE Ih , thus obtaining a twoby-two block system as in the middle form of Eq. 12. We then slightly alter the system by introducing the substitution E Ih = −E Ih and obtain the following discrete system of linear equations with the same definitions for the matrices as given in Eq. 11. In the sequel, we work with the system described in Eq. 13 unless otherwise stated, and equations and expressions involving matrices written in upper case bold letters denote problem-specific cases, whereas expressions with matrices written in calligraphic letters are valid in general.

Preconditioners for indefinite matrices of two-by-two block form
The matrix in equation (13) belongs to the class of indefinite matrices of two-by-two block form, where, in general A ∈ R n×n , B ∈ R m×n , C ∈ R m×n and D ∈ R m×m , n ≥ m. As pointed out in [23], under suitable partitioning, any linear system can be written in the twoby-two block-structured form (14) or in some equivalent form, such as For comprehensive studies of the properties of these matrices we refer to [23] and [7] and the references therein.
The experience how to precondition matrices of the form (14) is very rich. For the case when C = B we mention briefly some classical approaches, based on the exact block-LU factorisation of the matrix, namely, Here I n and I m denote identity matrices of size n and m, respectively. In, for instance, [23] and [13], the matrices very good preconditioners to the matrix A in equation (14), provided that we can construct A and S to be very good approximations of A and S and also accurately solve systems with those when applying the preconditioner. Matrices of the type (14) with square blocks occur not only when solving complex linear systems in real arithmetic, but also in various important applications, one example being the discrete first order necessary conditions in optimal control problems, constrained by partial differential equations. Below, we discuss preconditioners that utilise the same dimension of the diagonal and off-diagonal blocks.

Block-diagonal preconditioning
In a series of works (cf. [30,44,45,55,69,96,97]), based on some special norm techniques, another efficient preconditioner of block-diagonal form for matrices with square blocks has been derived and used. This preconditioner reads as follow for the general block systems If A and the symmetric part of B are symmetric and positive semi-definite matrices, the intersection of their nullspaces is empty and A + B is invertible. Furthermore, it has been shown that the spectrum of [79]). As a consequence, the condition number of the preconditioned system satisfies κ((P −1 bd 1 A) 2 ) ≤ 2 (cf. [10,30,96]). The above theoretical estimates show that this preconditioner is robust with respect to the discretisation and material parameters.
Applied to system Eq. 13, the preconditioning matrix is as follows For comparison purposes, in Section 4 we present numerical illustrations of the performance of P bd 1 for the target problem.

Schur complement-based preconditioning
In [63,64] for systems as in equation (13), efficient preconditioners have been derived, based on a particular high quality approximation ( S) of the arising exact (negative) Schur complement . The corresponding spectral bound for the Schur complement is shown to be Thus, the block-diagonal preconditioner reads as follows (cf. [64]), We see that the computation cost to solve a system with P bd 2 includes two solutions with M σ + (K − M ε ) and one solution with M σ , which is diagonal in the spectral element approximation for orthogonal hexahedral elements. Alternatively, one can use a block lower-triangular preconditioner with the same Schur complement approximation, The cost to use P bt 2 , compared to that when using P bd 1 is higher and requires one more matrix-vector multiplication with K − M ε . As noted in [10], this preconditioner has been shown to require slightly larger computational time compared to the PRESB preconditioner, outlined in Section 3.4 and therefore has not been tested in this study.

The PRESB preconditioner
Next, we consider the following preconditioner, referred to as 'PRESB', which stands for 'PREconditioning for Square Blocks' (cf. e.g., [10,14,15]). PRESB is derived for twoby-two block matrices (with square blocks) of the general form with the assumption that A is symmetric positive definite and that a and b are non-zero scalars having the same sign.
The theory covers also more general cases, for instance, when A is semi-definite, however these fall out of the scope of this paper. The preconditioner is of the form and it is shown (e.g., [10,14,15]) that all the eigenvalues of the preconditioned matrix P −1 A are real and are contained in the interval [ 1 2 , 1], independently of the discretisation parameter h, the parameters a and b and any other problem parameters that may be included in the matrix blocks. For completeness, we include the derivation of this result in Theorem 3.1 for a = b = 1 and more general types of matrices.
are all real and contained in the interval [ 1 2 , 1].
Since A and B are symmetric positive semi-definite and A + B is nonsingular, then all μ are real, nonnegative and less than 1. To estimate λ we use the following observation (with I being the identity matrix of appropriate dimension), From the latter expression we conclude that λ = 1 2 (1 + (1 − 2μ) 2 ). Using the bounds for μ, we obtain the following parameter-independent bounds for λ, It can easily be verified that inserting 0 and 1 for μ yields the upper bound for λ, whereas setting μ = 1 2 gives the lower bound for λ.
One additional advantage of the PRESB preconditioner is that it possesses the following block factorisation, For the problem, considered in this paper, we have that B 2 = B T 1 and a = b = 1. The preconditioner to A reads then From equation (23), we directly see the computational cost of applying the PRESB preconditioner, namely, solutions of linear systems with A + B and A + B T , one multiplication with the matrix B and three vector additions. For the matrix P the off-diagonal blocks are symmetric, and to solve a system of the form we solve twice with H = M σ +(K−M ε ), multiply with M σ and form three vector updates [10][11][12]. The computational procedure is shown in Algorithm 1.
The matrix H is of the size of the original complex matrix C A and can be also large. Therefore, in Steps 1 and 3 in Algorithm 1 we should also use an (inner) preconditioned iterative method. Thus, the overall solution procedure becomes of the so-called 'inner-outer' solution type. Due to the usage of an inner iterative solver the Algorithm 1 Solving linear system with preconditioner P.
PRESB preconditioner is of variable type and this imposes the additional requirement that the outer iterative Krylov subspace solver should be capable of handling variable preconditioning. Iterative methods that can be used as outer methods with variable inner preconditioning are the Flexible GMRES (FGMRES, [75]), the Generalized Conjugate Residual (GCR) method [34] and the Generalized Conjugate Gradient method [16]. In this study we use CGR as an outer solver and either the preconditioned GCR or FGM-RES methods for the inner solver with a suitable Algebraic Multigrid (AMG) preconditioning.
For clarity, the implementation of the iterative framework is summarised in Algorithm 2. In particular, the outer iterative solution procedure is given and solving the two inner systems in Algorithm 1 using an (inner) preconditioned iterative method refers to the inner solution method. Note that the outer loop is restarted every m + 1 iterations and that the current iterate x m+1 is taken as the new starting guess x 0 . It follows that at most m direction vectors p and m matrix-vector products with p need to be stored.
As can be seen, each outer iteration step involves solving a linear system with preconditioner P as described in Algorithm 1. Note that the right-hand side vector f used to solve system Eq. 25 corresponds to the normalised residual vector r of the current iteration of the outer loop.

Auxiliary-space Maxwell preconditioner
From Algorithm 1 and the form of the matrix blocks, it is seen that we have to solve systems with H = M σ + (K − M ε ). The matrix H is of type of a discrete form of the following variational problem arising from Maxwell's equations where c and d are scalars describing the magnetic and electric properties of the media, and where H 0,h is a subspace to H 0 (see Section 2.2 for the definition). We choose to solve linear systems of equations involving H by another (inner) preconditioned iterative solver. As a preconditioner we use a suitable AMG method. Efficient multigrid-based solvers for this problem are already available, described in detail in [48] and [53,54], and implemented in the library hypre [37,38] as the auxiliary-space Maxwell solver AMS [45,54]. The basic premise of this approach is based on the Hitpmair-Xu (HX) decomposition [48]. Without going into details, the underlying idea of their method lies in the understanding that the space H 0 (curl, ) possesses a stable decomposition (see Section 3 in [48] for a detailed discussion on decomposition of spaces) of form which implies that solving H 1 0 ( )-elliptic variational problems can be the basis for a preconditioner for H 0 (curl, )elliptic problems [48]. In essence, the auxiliary-space preconditioning technique comes down to approximately inverting the discrete curl-curl operator by transferring the problem into subspaces. Solving the problem in these (nodal) auxiliary spaces is much simpler than in the initial space as one can exploit the power of the AMG methods. As the auxiliary-space preconditioner can be viewed as a black box, we refrain from going into details and refer to [44,48,53,54] and [45] for details on the underlying theory and its implementation.
For completeness, we briefly specify the required user inputs for the auxiliary-space preconditioner. Besides the system matrix and right-hand side vector described by Eq. 26, AMS requires the commonly called discrete gradient matrix which describes the edges of the mesh in terms of its vertices [53,54]. If follows that the number of rows of this matrix corresponds to the number of degrees of freedom (edges) and the number of columns to that of the number of vertices in the mesh. In case of first order curl conforming basis functions, each row possesses two nonzero entries in the corresponding columns of the vertices composing the edge. In particular, the entries will be +1 and −1, where the sign is based on the orientation of the edge. In addition, the auxiliary-space preconditioner needs the coordinates of the vertices in the mesh to construct the solver. Alternatively, the user may also choose to input the representations of the constant vector fields in the edge element basis [53,54].
Note that the implementation of the auxiliary-space preconditioner AMS supports arbitrary order Nedelec elements as well as non-conforming quadrilateral/hexahedral meshes (i.e. meshes with hanging nodes).
We mention that the AMS implementation in hypre [37,38] holds for the positive semi-definite case [54] while convergence in the indefinite case might be slow or not reached at all. This could potentially pose some difficulties as, depending on the frequency as well as the air space conductivity considered, the system we solve might change type and become indefinite. For example, if the conductivity of air is set to 10 −8 S/m and the corresponding dielectric permittivity of air is 8.854 · 10 −12 F/m, then the system becomes indefinite in the air space of the domain at frequencies > 180 Hz due to the negative shift by M ε and the larger kernel of K [25]. Thus, for large frequencies the system matrix becomes highly indefinite.
In the numerical experiments, we use AMS as implemented in hypre [53,54] as the preconditioner for the iterative method applied to solve the inner systems (Steps 1 and 3 in Algorithm 1).

Numerical experiments
We illustrate the performance and the robustness of the solution procedure and the PRESB preconditioner for the EM equations on the following two test problems. Fig. 1 Table 1. Fig. 1(b). The conductive body which is inclined and of dimensions 1000× 662 × 3000 m 3 has an electric conductivity of 1 S/m. The surrounding host rock and the top layer have conductivities of 10 −4 and 0.01 S/m, respectively. The body is located at 2000 m depth and extends downwards to 5000 m depth. The top and the bottom of the structure stretch from −3000 to −4000 m and −2000 to −3000 m in x-direction respectively. In y-direction, the body extends from −331 to 331 m. The source is a grounded cable with its endpoints at coordinates (−75, 0, 0) and (59, 0, 0) m. The source moment is 100 Am. The computational domain is meshed with unstructured hexahedral elements (see Fig. 2) using the open-source meshing software gmsh [41]. Other relevant information pertaining to this setup is given in Table 1.

Problem 2 This problem is a 3D model that includes a conductive feature representative of, for example, an ore body buried within a resistive background and is covered by a thin conductive layer as shown in
The numerical procedures described in Sections 2 and 3 are implemented in a parallel fashion using the Message Passing Interface (MPI) and employ the open-source libraries PETSc [20,21] and hypre [37,38].
All experiments are performed on AMD Ryzen Threadripper 2950X 16-core processor with a clock frequency of 3.5 GHz and with 128 GB RAM. The operating system is Ubuntu 20.04. All numerical simulations are run with two MPI processes in order to show that the developed iterative framework and all libraries used within it can run in parallel fashion.
As described in Section 3, when the system is solved iteratively, the outer iteration method for the two-by-two block system is the GCR method. The inner two systems arising from applying the preconditioner are solved using either GCR or FGMRES provided by PETSc. In both cases, AMS is used as a preconditioner. For comparison we provide experiments, where the inner systems are solved by the sparse direct solver MUMPS, also available through PETSc.
We remark that we also tested hypre's ILU called pilut as a potential preconditioner for the inner systems. However, we do not include any results of these simulations as all of them failed to reach convergence for Problems 1 and 2. We thus note that ILU decomposition does not constitute a suitable inner preconditioner for the system to be solved in the proposed iterative framework.
For all presented tests, the electric conductivity of the air is set to 10 −8 S/m. The magnetic permeability and dielectric permittivity are assumed to be the corresponding constants of free space for all but a few indicated exceptions. Unless stated otherwise, by default the simulations are run using the smallest problem size (see Table 1) and are stopped once the relative residual of the outer solver GCR falls below the threshold of 10 −12 . Regarding the stopping tolerance for the inner GCR solver preconditioned with AMS, experiments show that it can be chosen to be  Table 2 indicates that the stopping tolerance for the inner solver can be set to 10 −3 while still allowing for a good compromise between the number of outer and inner iterations to keep the solver robust and the computational time down. Further, we note that the given numbers of degrees of freedom (DOF) correspond to the real-valued system, which is twice that of the complexvalued system.
The focus of this work is on the convergence behaviour of the described iterative framework. Nonetheless in order to demonstrate that the iterative algorithm has the required modelling accuracy compared to other established methods, a numerical comparison of the electromagnetic responses to finite element solutions for the 3D model (Problem 2) is shown in Fig. 5 in Appendix A. Details are disclosed in Appendix A.
The conductivity distributions of Problems 1 and 2 shown in Fig. 1 involve several difficulties that might negatively affect the convergence of the used iterative methods, such as large conductivity contrasts, in particular at the air-Earth interface and to a lesser extent within the Earth at interfaces between different entities, non-uniformly stretched cells, unstructured elements and locally refined meshes. Hence, on the basis of these problems, we ascertain the robustness of the solvers.
As predicted by the convergence properties of the PRESB preconditioner, we observe that across the tested frequency range of 0.1 to 10'000 Hz the iteration counts of the outer solver remain virtually stable (see Tables 3 and 4 as well as Fig. 3). It is worth noting that solving the two inner systems (Steps 1 and 3 in Algorithm 1) iteratively leads to only a slight increase of the outer iteration count compared to the corresponding count where the inner systems are solved by a direct method and thus exactly (compare outer iteration counts in Tables 3 and 4). Another observation is that solving the two inner systems iteratively requires considerably less memory than when employing a direct solver for the same two inner systems. It is also evident that using preconditioned iterative solution techniques for solving the two inner systems instead of a direct solver reduces the computational time spent to reach convergence with the outer iterative algorithm. This gain in time only manifests if the preconditioned iterative algorithm applied to the two inner systems is efficient which in fact means the applied solution technique needs to converge to the desired tolerance in few iterations. In Tables 3 and 4, we observe that the inner preconditioned iterative method performs effectively for all but the highest (10'000 Hz) tested frequencies. In this case, the solution time increases drastically (approximately by a factor of five compared to the next lower frequency of 8000 Hz) while the outer iteration count remains unaffected which points to a break down of the efficiency of the inner preconditioned iterative solver. This is indeed the case as visualised in Fig. 4, which shows the convergence of the inner preconditioned iterative solver for the two inner systems at two stages where the residual of the outer algorithm dips below 10 −4 and 10 −8 , respectively. It can be seen that for a frequency of 10'000 Hz the inner solver needs significantly more iterations to reach  convergence than for the other frequencies, which explains the observed increase in overall solution time. Thus, to achieve good overall performance we need a good solver for the inner systems. As a side note, clearly, the deterioration of the performance of the inner iterative solver can be detected and, if memory consumption allows to do so, one might switch to solving the inner systems directly. As alluded to in Section 3.4.1, the inner system H becomes indefinite for frequencies larger than 180 Hz and, as already noted, this is due to the negative shift from M ε and the large nullspace of the stiffness matrix K. This affects negatively the performance of the preconditioner for the inner solver as AMS is not devised for indefinite systems (cf. [54]). This effect can clearly be observed for the highest frequency of 10'000 Hz in Fig. 4.
In this study, we use AMS-preconditioned GCR as an inner solver and we illustrate that the AMS preconditioner is not efficient enough as soon as the inner systems become strongly indefinite. However, we see that, as theory predicts, the outer convergence rate does not get affected because of that. Thus, a viable option for higher frequencies is to use another inner solver, that is more suitable for such inner systems.
The experiments convincingly show that the convergence of the outer solver is independent of the problem sizethe number of iterations of the outer solver remains nearly constant and the solution time increases linearly with increasing number of degrees of freedom, see Table 5. For additional insight in the performance of the PRESB preconditioner, in Table 6 we compare the solution times and peak memory requirements, when the inner systems are solved both iteratively and by a direct method and those when solving the original complex-valued system in Eq. 10. As expected, the iterative solver has significantly lower memory demands and is faster than when direct methods are used. For simulations in which either the complex system or  the inner systems are solved directly one can easily spot the increase in memory consumption as well as in solution time.
The execution time given in the tables include the overall time for the iterative solver to converge, including the time needed to construct the AMS preconditioner. Whenever direct methods are used, the analysis and factorisation times are included. For complex systems we report the time for analysis, factorisation and solution time. Memory requirements in particular may become the limiting factor for large problems. This is observed for the largest problem size where the direct solvers demand approximately 151 and 198 GB of memory respectively, which exceeds the available memory of the computer platform used to perform these tests. In contrast, using an iterative method for the inner systems uses only around 30 GB of memory. In order to test the robustness of the solver with respect to magnetic permeability μ, we assign different values of relative magnetic permeability μ r to the ore body in Problem 2. In particular, μ r is set to 2, 5 and 10. Per definition, the relative magnetic permeability is defined as μ r = μ μ 0 . We point out that the magnetic permeability for different rock types falls within a very narrow range between 1μ 0 and 3μ 0 . Further, a relative magnetic permeability value of 10 can safely be considered to cover the vast majority of possibilities encountered in geo-EM applications. The test results are presented in Table 7 and show that the influence of spatially variable magnetic permeability distributions is minimal thus confirming that the solver is robust also with respect to material parameter μ.
We next consider the robustness of the solver with regard to variable dielectric permittivity ε. To do so, we assume that the surrounding host rock around the buried body and the cover layer on top of the host rock have relative dielectric permittivities values of 5 and 20, respectively. In a geophysical sense, this could be representative of a weakly fracture granite covered by some kind of glacial deposits. Further, the relative magnetic permeability of the ore body is taken to be either 1 or 10. The results for this set of parameter value combinations are given in Table 8, and, in comparison to Table 7, these imply that the solver remains robust with regard to variable dielectric permittivity distributions.  Taken all together, our examples demonstrate that the solver is robust with respect to all material properties σ , μ, and ε relevant in geophysical settings and frequency ω as well as with respect to the discretisation parameter. Moreover, the outer solver convergences within relatively few iterations. The whole iterative framework is computationally inexpensive in terms of memory requirement when compared to direct solvers. On top of that, the procedure also saves a considerable amount of computational time given that the inner systems can be solved efficiently (i.e. in a few iterations).

Convergence behaviour for loop sources
In this section, we test the convergence behaviour of our algorithm for Problem 1 when using a horizontal or vertical loop source corresponding to vertical or horizontal magnetic dipole, respectively, instead of line sources. The loop sources are square loops with a side length of 10 m, and the centres of the horizontal and vertical loops are located at coordinates (0, 0, 0) and (0, 0, 5), respectively. Each side is subdivided into four wire segments. We remark that the simulations for the loop sources are performed using the intermediate problem size with 3'641'400 degrees of freedom. Tables 9 and 10 show the relevant convergence parameters, namely the outer iteration count and the time taken to reach convergence when solving all the systems in the algorithm iteratively. For comparison, we also run the simulations where the inner systems are solved using a direct solver. As already pointed out earlier, solving the inner systems by iterative methods is more cost-effective, both in terms of time as well as memory required to carry out the simulation. However, we observe that the more complicated nature of the loop source that enters the problem through the right-hand side increases the overall simulation times when compared to those of a line source (see Table 5). This is due to an increase in iterations required to reach convergence for the inner systems as demonstrated by the average number of inner iterations per system given in Tables 9 and  10. We recall that the increase in time can be compensated for by relaxing the inner stopping criteria which results in fewer inner iterations at the expense of more outer iterations (see Table 2). In the case of loop sources, this might become a viable option to keep the simulation times low. Even though the time gain of solving the inner systems by a preconditioned iterative method is not as large as for the same simulations with line sources, the overall memory consumption still presents a major advantage over using direct solvers as the method of choice to solve the inner systems.

Comparisons of the performance of P bd 1 and P
We next present a comparison of the PRESB preconditioner P with the block diagonal preconditioner P bd 1 , used in related works (e.g. [30,44,45,69]). For the sake of clarity, we briefly outline the general setup of the solver. Recall the two-by-two system to be solved preconditioned with a block diagonal matrix of the form We note that none of the aforementioned works (namely [30,44,45,69]) include displacement currents, that is M ε is present neither in their block systems nor in the used    block preconditioners. In addition, the outer solver used in these publications is either a GMRES or FGMRES, and the two inner systems are solved using a AMS-preconditioned CG. Hence, the choices of the outer and inner solvers differ, however GCR and GMRES are shown to be mathematically equivalent, (cf. e.g. [50]). Thus, we can make a fair comparison by testing the performance of the block diagonal preconditioner for our problems. Table 11 shows the performance of both block preconditioners in terms of outer iteration count (N outer it ), time spent inside the algorithm (time [s]) and peak memory consumption (mem[GB]) for each simulation across a frequency range spanning five orders of magnitude. It is easily discernible that both preconditioners exhibit similar behaviour: (1) an increase of outer iterations for higher frequencies; (2) highest iteration count for the frequency of 100 Hz; (3) a steep rise in time spent to reach the solution for frequency 10'000.
The comparison clearly shows that the PRESB preconditioner is more efficient than the block-diagonal preconditioner, even though, compared to the application cost of the block-diagonal preconditioner, applying it requires one additional matrix-vector multiplication and four vector update operations. Note, in particular, that the outer iteration counts when using the PRESB preconditioner are equal to or lower than when applying the block-diagonal one. The same behaviour also holds true (except for frequency 100 Hz) when reducing the outer stopping tolerance to 10 −8 for simulations in which the block-diagonal preconditioner is used. As the outer iteration count is smaller for runs using PRESB, the time spent to achieve convergence is lower.
All in all, given the presented evidence, the PRESB preconditioner is to be preferred. As the implementation of both preconditioners is fairly similar, changing to the improved preconditioner is beneficial and thus encouraged.

Conclusions and outlook
In this work, we present an efficient and robust iterative framework for solving the linear algebraic system arising from the spectral element discretisation of the curl-curl equation of the total electric field describing the behaviour of the electromagnetic field for 3D CSEM problems in frequency domain. We transform the original complexvalued system into an equivalent real-valued form. Then,  we solve the resulting real-valued two-by-two block system by the GCR method preconditioned with a highly efficient block-based preconditioner PRESB. Our numerical experiments confirm the theoretical estimates that the PRESB-preconditioned GCR solver is robust with respect to spatially variable material parameters σ , μ and ε and with respect to the number of degrees of freedom by showing experiments across five orders of magnitude with regard to problem size and for stretched structured, unstructured as well as locally refined meshes. For the models considered in this work, the iterative algorithm reaches convergence within a small (typically < 20) number of outer iterations. The efficiency of the iterative framework depends on the efficiency of the solution of the two inner systems of equations, being a part of PRESB. In this study, these two inner systems are either solved using a direct method or an iterative GCR/FGMRES method with the auxiliaryspace preconditioner AMS. Our numerical tests clearly show that it is computationally beneficial to employ the suggested preconditioned iterative solver to the two inner systems instead of a direct solver in terms of overall memory requirements and time spent to obtain the solution for the two-by-two block system. Further and even though the iterative framework works with a system twice the size of the original complex-valued system, the iterative method works on large-scale problems where direct solvers succumb due to their excessive memory consumption, even on smaller less performance oriented computers and platforms.
From the above, it becomes clear that the efficiency of the overall algorithm stands and falls with the performance of the auxiliary-space preconditioner which is devised for semi-definite H 0 (curl, ) problems. We note that the system matrix of the two inner systems does not necessarily adhere to semi-definiteness. In fact, the system matrix becomes indefinite for frequencies larger than 180 Hz. Our experiments show that the performance of AMS loses its advantages in terms of time over the direct solver for a frequency of 10'000 Hz. Thus, to achieve the full efficiency of PRESB across all possible frequency ranges, we need to use either frequency-dependent inner solvers or an advanced variant of AMS that could ensure better efficiency in the indefinite case.
Comparison to a similar iterative framework preconditioned with a block-diagonal preconditioner reveals that our approach converges in slightly fewer iterations resulting in certain execution time gains. It is noted that previous works show strong scalability for the iterative solver as well as for AMS. Further, AMS provides support for high-order curlconforming element discretisations. Our future research is thus focused on extending our framework for high-order basis functions and to upscale our code in order to exploit the performance offered by modern distributed-memory platforms.

Appendix A: Accuracy for the 3D model compared to finite element solutions
In the following, we present results from numerical simulations obtained by our iterative framework and by a finite element code developed by [74], that uses a direct solver. In particular, we deploy twelve receivers across the target body of the 3D model (Problem 2) in cross-line configuration, meaning the receiver line is perpendicular to the transmission direction. The receiver line is offset by by 3.5 km in x-direction. The receivers range from −3 to −0.5 km and from 0.5 to 3 km in y-direction at intervals of 500 m. Figure 5 shows the electric and magnetic fields excited by line source transmitting at 10 Hz for both the iterative framework and the finite element code denoted by SEM and FEM respectively. The responses for both numerical schemes are in good agreement with each other. Further, Table 12   and magnetic field components across all receivers thus confirming the observed agreement. Note, that the relative deviations are given for frequencies of 10 and 100 Hz. 2018-05973.
Funding Open access funding provided by Uppsala University.

Data Availability
The data underlying this article can be shared upon reasonable request. Since the modelling code is to be embedded in a larger software package, it cannot be made freely available.

Conflict of Interests
The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.