With application to acoustic, elastic and electromagnetic scattering, the geological formations often can be approximated by horizontally layered background media with several embedded inhomogeneities. For 3D case, such problems usually lead to the large ill-conditioned system of linear equations when finite differences (or finite elements) are employed for discretization. In this work we assume that the discretization grid
is a tensor product of one-dimensional grids. Direct solvers can be efficiently applied for the solution of 2D (or 2.5D) problems as well as 3D problems of rather small size only. Therefore, the iterative solvers are the only option for 3D large-scale problems. The choice of preconditioner is crucial for robustness of iterative method. In [8, 9] we show the efficiency of the preconditioning via discrete Green’s function in layered medium. Similar ideas are used in the volume integral equation (IE) for the scattered field . However, a well known drawback of the IE is the large cost of numerical integration. In particular, provided the horizontal layering and uniform discretization grid in and , FFT is applied along and and correlation is computed along , so the cost is . Here , and are numbers of discretization nodes along , and , respectively. For non-uniform discretization in and , discrete Fourier transform is employed and the cost is .
Quadratic dependence of the cost on was removed in FDIE approach of [8, 9]. FDIE can be viewed as a finite-difference discretization of the volume IE approach. But thanks to the sparsity of the finite-difference discretization, the cost of preconditioner in FDIE is and for uniform and non-uniform grids, respectively (similar approach was proposed in  for volume IE method).
For geophysical applications one typically uses fine (not necessarily uniform) grid within the domain of the survey and coarse (exponentially convergent) grid outside (possibly, with imaginary grid steps in PML). Therefore the cost of preconditioners of [8, 9, 1] is still quadratic with respect to and . In this paper we propose a fast algorithm to compute the solution of Helmholtz equation in layered medium for arbitrary right hand side part. The algorithm is based on the ideas of cyclic reduction (see [3, 4]) which is is closely related to specific method of domain decomposition type. We split the computational grid into connected Cartesian blocks and compute approximate solution by solving the problem on each block independently with homogeneous Dirichlet boundary conditions. The corresponding residual has support on the interfaces of the blocks . Therefore, by performing just local (on ) computations we obtain Dirichlet conditions for on each interface . The last step of our algorithm is computing in each under given boundary conditions . The cost of the algorithm is and for uniform and non-uniform grids in the domain of the survey, respectively. We also extended our approach for solving Maxwell’s equation in the frequency domain. In a very similar way the algorithm can be applied for solving elasticity problems as well.
2 Scalar problem
We consider the Helmholtz equation in horizontally layered medium
in , where , is regular enough positive function of .
Let , and be (perhaps, complex if PML is considered for ) grid steps of , and , respectively. Consider 7-point finite-volume discretization scheme
where are lengths of the edges of the control volume .
The system (2) can be rewritten as
where , and are discretizations of , and , respectively, and , and are diagonal mass matrices: , for example.
Split the computational grid into connected Cartesian blocks and define and the projection operators to interior nodes of and , respectively. Let and be projected operators and , where . In each block we construct the approximate solution satisfying the equation
and homogeneous Dirichlet boundary conditions . Here . We also define the global approximate solution as and . Denote and:
First we Fourier-transform (2) with respect to and :
where , . Obviously, the system (2) can be split (with respect to eigenmode) into multiple tridiagonal systems. For each pair and each eigenmode the solution can be computed in operations.
It is easy to see that the residual has a support only on . Therefore, thanks to the sparsity of , to obtain we need only on as well as at neighboring nodes. Under known , this step can be performed via inverse Fourier transform.
Obviously, the difference satisfies the homogeneous equation in each with some Dirichlet boundary conditions . Thanks to the sparsity of , the latter can be obtained in operations via Fourier transform of (2). In fact, let and be diagonal matrix of generalized eigenvalues and matrix of generalized eigenvectors of operator , respectively: . Then satisfies
where Fourier transform can be performed in operations thanks to the sparsity of . The system (2) can be split into tridiagonal systems and each of them can be solved in operations. Finally, to obtain we need to perform another operations of inverse Fourier transform on the .
To compute in the interior nodes of each , we first solve the Fourier transformed equation
in , where satisfies the computed boundary conditions . Finally, in each we perform inverse Fourier transform and obtain .
Our algorithm can be summarized as follows
for each perform Fourier transform and compute
for each solve the problem (2) and compute
via inverse Fourier transform, compute at the nodes adjacent to .
compute residual at
compute via Fourier transform with respect to and
under known , compute for each
compute via inverse Fourier transform of in each
The computational costs of the second, the fourth and the fifth steps are , and , respectively. The costs of the remaining steps of the algorithm depend on how we split our grid into blocks . For the case of uniform grids in the domain of the survey we split each of into three subgrids: and correspond to the non-uniform exponentially convergent grids and covers the uniform part. To maintain the same approximation error within the domain of the survey and outside, further we assume that the number of nodes in exponentially convergent grids along and is proportional to logarithm of uniform parts of and , respectively. For the case of general non-uniform grids we split () into subgrids with nodes in each of them. For the first type of gridding we compute Fourier transforms in uniform parts via FFT. Therefore, the costs of the first, the sixth and the seventh steps are . For general non-uniform gridding, the cost of Fourier transforms on each is and the total number of blocks is . Hence, the overall costs of the first, the sixth and the seventh steps are . The cost of the third step is and for two type of griddings above.
By extending our approach to multi-level cyclic reduction (i.e. by employing multiple mutually embedded partitions of ), we can reduce the asymptotics of the cost from to along each direction. But the constant factor is significantly greater due to multiple (on each level of embedding) solution of block-tridiagonal systems.
3 Maxwell equations
Consider the Maxwell equation for magnetic field in the layered medium:
Here we made an assumption of negligible displacement current which is typical for large-scale problems of geophysical prospecting.
We consider a bounded computational domain and introduce a Cartesian three-dimensional (3D) grid as follows:
where , , and are the number of grid nodes in the , and directions. The Lebedev -grid is defined as a sub-grid of with even sum of indices and the Lebedev -grid contains the remaining nodes of . For example, for even the grid in the plane is shown in Fig. 1.
In other words, each of the Lebedev - and -grids is a combination of two Yee grids in the two-dimensional (2D) configuration and four Yee grids in 3D configuration (see Fig. 1). For any one-dimensional (1D) grid we define primary and dual grids as sub-grids of
with even and odd indices, respectively. Therefore, in theplane and for even k, the Lebedev -grid is a superposition of the tensor product of two primary grids as well as the tensor product of two dual grids. On the other hand, for odd , the Lebedev -grid is a superposition of the tensor product of a primary (along ) and a dual (along ) grid as well as the tensor product of a dual (along ) and a primary (along ) grid. In fact, according to the definition of the -grid and for even , it consists of vertices with either both even and (which is a tensor product of two primary grids) or both odd and (which is a tensor product of two dual grids). Similarly, for odd we have a tensor product of primary and dual grids as well as a tensor product of dual and primary grids.
Let all three components of the finite-difference vector magnetic fieldbe defined at the same nodes of the -grid, and similarly, all three components of the finite-difference vector electric field be defined at the same nodes of the
-grid. Hence there is no need to perform any interpolation with our scheme to handle anisotropic Ohm’s law. We define finite-difference derivatives on both grids along thedirection as follows:
and similarly along and direction. Note that equation (9) performs mapping from -grid to -grid and vice versa. This allows us to obtain a self-consistent discrete system of Maxwell equations
here we use the same notations for discrete fields as for continuous ones.
After eliminating from (3) we obtain the linear system of equations for the magnetic field vector
The boundary of the computational domain contains nodes of both - and -grids. We replace the boundary conditions at infinity by:
As pointed out in , this allows us to decrease the computational domain by a factor of 2 (compared to the conventional boundary conditions or ) without losing accuracy.
Consider the Fourier transform (with respect to and ) of equations (3). Due to presence of mixed derivatives in (3), the equations in Fourier domain can not be represented in such a simple Krononecker product form as (2). According to the definition of the Lebedev grid, for even grid node in the -plane, the magnetic field is defined on a tensor product of two primary grids and a tensor product of two dual grids. Denote them by and , where the first and second indices mean the type of the grid along and , respectively. Similarly, and are defined for odd
. Define eigenvalues and eigenfunctions of 1D operators on each grid with Dirichlet conditions on primary grid and Neumann conditions on dual grid:
where . One can derive the relations between primary and dual eigenfunctions: and . For each , we expand and as a sum
Here for even and for odd . In this case, the equations for the expansion coefficients have the following form (we omit the indices and of , , and to be more concise):
For even and for the primary-primary grid:
Here and below and finite-difference derivatives are defined by (9). For even and for the dual-dual grid:
For odd and for the primary-dual grid:
Finally for odd and for the dual-primary grid: