A nested Schur complement solver with mesh-independent convergence for the time domain photonics modeling

TL;DR

This paper introduces a nested Schur complement solver for Maxwell equations that achieves mesh-independent convergence, improving the efficiency of time domain photonics simulations involving complex boundary conditions.

Contribution

It presents a novel nested Schur complement approach that effectively solves double saddle point systems with mesh-independent convergence in photonics modeling.

Findings

Achieves mesh-independent convergence in the outer solver

Efficiently handles elliptic inner systems with various solvers

Applicable to exponential and implicit time integration methods

Abstract

A nested Schur complement solver is proposed for iterative solution of linear systems arising in exponential and implicit time integration of the Maxwell equations with perfectly matched layer (PML) nonreflecting boundary conditions. These linear systems are the so-called double saddle point systems whose structure is handled by the Schur complement solver in a nested, two-level fashion. The solver is demonstrated to have a mesh-independent convergence at the outer level, whereas the inner level system is of elliptic type and thus can be treated efficiently by a variety of solvers.

Tables

Table 1. Table 1: The number of degrees of freedom for the meshes used in the tests

mesh	system size
$n_x\times n_y\times n_z$	$N=n+m$
$20\times 20\times 12$	$\mathrm{45\hspace{0.17em}565}=\mathrm{34\hspace{0.17em}398}+\mathrm{11\hspace{0.17em}167}$
$40\times 40\times 24$	$\mathrm{333\hspace{0.17em}425}=\mathrm{252\hspace{0.17em}150}+\mathrm{81\hspace{0.17em}275}$
$80\times 80\times 48$	$\mathrm{2\hspace{0.17em}548\hspace{0.17em}441}=\mathrm{1\hspace{0.17em}928\hspace{0.17em}934}+\mathrm{619\hspace{0.17em}507}$
$160\times 160\times 96$	$\mathrm{19\hspace{0.17em}922\hspace{0.17em}345}=\mathrm{15\hspace{0.17em}086\hspace{0.17em}022}+\mathrm{4\hspace{0.17em}836\hspace{0.17em}323}$

Table 2. Table 3: Results for the nested Schur complement solver for the “difficult” test case, with random exact solution vector.

mesh	residual	iterations	CPU
	norm	outer (inner)	time
$20\times 20\times 12$	7.50e-11	31 (68)	3.59 s
$40\times 40\times 24$	6.67e-11	32 (108)	17.3 s
$80\times 80\times 48$	6.84e-11	32 (145)	135.6 s
$160\times 160\times 96$	7.01e-11	31 (234)	1258 s

Table 3. Table 4: Comparison results of the nested Schur complement solver and the FS preconditioner. The former is employed with GMRES(10) and ICCG(0) as the outer and inner solvers, respectively. The FS preconditioner is applied with nonrestarted GMRES.

mesh,	method	CPU time,	residual
tolerance		iter outer(inner)	norm
$40\times 40\times 24$	FS prec.	0.49 s, 7	7.58e-12
9.64e-11	nested Schur	0.54 s, 2(8)	1.23e-13
$80\times 80\times 48$	FS prec.	3.84 s, 8	5.95e-10
8.09e-09	nested Schur	3.51 s, 2(8)	2.09e-09
$160\times 160\times 96$	FS prec.	59.7 s, 14	3.97e-09
4.23e-09	nested Schur	45.8 s, 2(19)	3.77e-09

Equations

$$\{\begin{array}{rlrlrlrl}& \mu {\partial }_{t}H=& -{\sigma }_1& H& -\nabla \times & E& +& J_H,\\ & \epsilon {\partial }_{t}E=& \nabla \times & H& -{\sigma }_2& E& +& J_E,\end{array}$$

$$\left[\begin{array}{cc}{M}_{\mu }& 0\\ 0& M_{\epsilon }\end{array}\right]\left[\begin{array}{c}{\mathbf{h}}^{\prime }\\ {\mathbf{e}}^{\prime }\end{array}\right]=-\left[\begin{array}{cc}{M}_{{\sigma }_1}& K\\ -K^T& M_{{\sigma }_2}\end{array}\right]\left[\begin{array}{c}\mathbf{h}\\ \mathbf{e}\end{array}\right]+\left[\begin{array}{c}{\mathbf{j}}_H\\ {\mathbf{j}}_E\end{array}\right],$$

$$\begin{array}{c}\left[\begin{array}{c}{\mathbf{h}}^{\prime }\\ {\mathbf{e}}^{\prime }\end{array}\right]=-A\cdot \left[\begin{array}{c}\mathbf{h}\\ \mathbf{e}\end{array}\right]+\left[\begin{array}{c}{M}_{\mu }^{-1}{\mathbf{j}}_H(t)\\ M_{\epsilon }^{-1}{\mathbf{j}}_E(t)\end{array}\right],\\ A=\left[\begin{array}{cc}{M}_{\mu }^{-1}& 0\\ 0& M_{\epsilon }^{-1}\end{array}\right]\left[\begin{array}{cc}{M}_{{\sigma }_1}& K\\ -K^T& M_{{\sigma }_2}\end{array}\right]=\left[\begin{array}{cc}{M}_1& K_1\\ -K_2^T& M_2\end{array}\right]\in {\mathbb{R}}^{n\times n},\end{array}$$

$$y^{\prime }(t)=-\mathcal{A}y(t)+g(t),\mathcal{A}=\left[\begin{array}{cc}A& B_1^T\\ -B_2& 0\end{array}\right]\in {\mathbb{R}}^{N\times N},$$

$$(I+\gamma \mathcal{A})x=b,$$

$$\left[\begin{array}{cc}\mathtt{A}& \mathtt{B}\\ \mathtt{C}& \mathtt{D}\end{array}\right],$$

$$\left[\begin{array}{cc}\mathtt{A}& \mathtt{O}\\ \mathtt{O}& \mathtt{D}-\mathtt{C}{\mathtt{A}}^{-1}\mathtt{B}\end{array}\right]\text{or}\left[\begin{array}{cc}\mathtt{A}-\mathtt{B}{\mathtt{D}}^{-1}\mathtt{C}& \mathtt{O}\\ \mathtt{O}& \mathtt{D}\end{array}\right].$$

$$\mathtt{D}-\mathtt{C}{\mathtt{A}}^{-1}\mathtt{B}=I+{\gamma }^2{B}_2(I+\gamma A{)}^{-1}{B}_1^T$$

$$\mathtt{A}-\mathtt{B}{\mathtt{D}}^{-1}\mathtt{C}=I+\gamma A+{\gamma }^2{B}_1^T{B}_2$$

$$I+\gamma A=\left[\begin{array}{cc}I+\gamma M_1& \gamma K_1\\ -\gamma K_2^T& I+\gamma M_2\end{array}\right]$$

$$I+\gamma M_2+{\gamma }^2{K}_2^T(I+\gamma M_1{)}^{-1}{K}_1,$$

$$I+\gamma M_1+{\gamma }^2{K}_1(I+\gamma M_2{)}^{-1}{K}_2^T.$$

$$I+\gamma M_2+{\gamma }^2{K}_2^T(I+\gamma M_1{)}^{-1}{K}_1=M_{\epsilon }^{-1}\left[M_{\epsilon }+\gamma M_{{\sigma }_2}+{\gamma }^2{K}^T(M_{\mu }+\gamma M_{{\sigma }_1}{)}^{-1}K\right],$$

$$\tilde{\mathcal{A}}\tilde{x}=b,\tilde{\mathcal{A}}=(I+\gamma \mathcal{A}){\mathcal{M}}^{-1},\tilde{x}=\mathcal{M}x,$$

$$\mathcal{M}=(I+\gamma {\mathcal{A}}_1)(I+\gamma {\mathcal{A}}_2),{\mathcal{A}}_1+{\mathcal{A}}_2=\mathcal{A}.$$

$$\begin{array}{c}{\mathcal{A}}_1=\left[\begin{array}{cc}{A}_1& B_{1,H}^T\\ -B_{2,H}& 0\end{array}\right],{\mathcal{A}}_2=\left[\begin{array}{cc}{A}_2& B_{1,E}^T\\ -B_{2,E}& 0\end{array}\right],\\ A=A_1+A_2,A_1=\left[\begin{array}{cc}{M}_1& K_1\\ 0& 0\end{array}\right],A_2=\left[\begin{array}{cc}0& 0\\ -K_2^T& M_2\end{array}\right],\end{array}$$

$$B_1=B_{1,H}+B_{1,E},B_2=B_{2,H}+B_{2,E},$$

$${\hat{B}}_1^T=\left[\begin{array}{cccc}{K}_1& 0& -I& 0\\ 0& -K_2^T& 0& -I\end{array}\right],{\hat{B}}_2=\left[\begin{array}{cc}0& \Sigma \\ \Sigma & 0\\ -{\Sigma }_{\ast }& 0\\ 0& -{\Sigma }_{\ast }\end{array}\right],$$

$$\Sigma =\mathrm{diag}\left[\begin{array}{ccc}{\mathbf{\sigma }}_x,& {\mathbf{\sigma }}_y,& {\mathbf{\sigma }}_z\end{array}\right],{\Sigma }_{\ast }=\mathrm{diag}\left[\begin{array}{ccc}{\mathbf{\sigma }}_y\cdot {\mathbf{\sigma }}_z,& {\mathbf{\sigma }}_x\cdot {\mathbf{\sigma }}_z,& {\mathbf{\sigma }}_x\cdot {\mathbf{\sigma }}_y\end{array}\right].$$

$${\hat{B}}_1^T{\hat{B}}_2=B_1^T{B}_2=\left[\begin{array}{cc}{\Sigma }_{\ast }& K_1\Sigma \\ -K_2^T\Sigma & {\Sigma }_{\ast }\end{array}\right],$$

$${\hat{B}}_1={\hat{B}}_{1,H}+{\hat{B}}_{1,E}{\hat{B}}_2={\hat{B}}_{2,H}+{\hat{B}}_{2,E},$$

$${\hat{B}}_{1,H}^T:=\left[\begin{array}{cccc}{K}_1& 0& -I& 0\\ 0& 0& 0& 0\end{array}\right],{\hat{B}}_{2,H}:=\left[\begin{array}{cc}0& \Sigma \\ 0& 0\\ -{\Sigma }_{\ast }& 0\\ 0& 0\end{array}\right],$$

Full text

1 Introduction

Numerical solution of the time-dependent Maxwell equations is an important computational problem arising in various scientific and engineering fields such as photonic crystal modeling, gas and oil industry, biomedical simulations, and astrophysics. Rather often the application environment suggests that the Maxwell equations have to be solved many times, for instance, for different source functions or different medium parameters) [6]. The size of the spatial computational domain as well as the necessity to solve the equations many times make this task very demanding in terms of computational costs. Therefore, advanced computational techniques have to be applied, such as modern finite element discretizations [9, 28] in space and efficient integration schemes in time. Along with multirate and implicit time integration schemes [31, 10], exponential time integration schemes [18] have recently been shown promising for solving the Maxwell equations in time [19, 5, 6].

Exponential time integration schemes, which are essentially based on the notion of the matrix exponential, are attractive not only because of their excellent stability and accuracy properties but also due to their efficiency and potential for a parallelism in time [15, 21]. Most frequently, especially when the spatially discretized Maxwell operator $\mathcal{A}$ (defined in (3) below) is not a skew-symmetric matrix, the actions of the matrix exponential in exponential time integration schemes are evaluated by Krylov subspace methods. To be efficient, Krylov subspace methods often need to rely on rational approximations [11, 23, 30, 16] (so that the Krylov subspace is built up for a rational function of $\mathcal{A}$ rather then for $\mathcal{A}$ itself) and on the so-called restarting techniques [25, 29, 1, 16, 12] (to keep the Krylov subspace dimension restricted). A popular variant of the rational Krylov subspace methods is the shift-and-invert (SAI) method [23, 30]. Rational Krylov subspace methods and, in particular, the SAI Krylov subspace method as well as implicit time integration schemes involve solution of linear systems with the matrix $I+\gamma\mathcal{A}$, with $\gamma>0$ being a given parameter (which is, in case of implicit time stepping, the time step size).

Despite the significant progress achieved last decennia in sparse direct linear system solvers, for three-dimensional (3D) problems iterative linear system solvers remain the methods of choice. The task of solving linear systems with the matrix $I+\gamma\mathcal{A}$ when $\mathcal{A}$ is a spatially discretized Maxwell operator is especially challenging for the Maxwell equations. This is caused not only by the fact that the matrix $\mathcal{A}$ has a saddle point structure but also due to special nonreflecting boundary conditions which are often imposed for the Maxwell equations. In this paper a preconditioner is proposed to solve iteratively linear systems with the matrix $\mathcal{A}$ when the popular perfectly matching layers (PML) boundary conditions are imposed. In this case the matrix $\mathcal{A}$ has the so-called double saddle point structure [2].

This paper is organized as follows. The problem is set up in Section 2. In Section 3 we present the nested Schur complement solver. Other possible preconditioners for these problems are discussed Section 4. Section 5 is devoted to numerical experiments. Finally, the conclusions are drawn in the last section whereas some background material is given in two appendices.

2 Problem formulation

We are interested in solving a system of time-dependent three-dimensional (3D) Maxwell equations

[TABLE]

where $H=H(x,y,z,t)$ and $E=E(x,y,z,t)$ are respectively magnetic and electric fields, $\mu$ is the magnetic permeability (as typical for photonics and gas-and-oil exploration applications, $\mu\equiv 1$ for all the tests considered in this paper; however, in general one can have $\mu=\mu(x,y,z)$) and $\varepsilon=\varepsilon(x,y,z)>0$ is the electric permittivity. Furthermore, $\sigma_{1,2}=\sigma_{1,2}(x,y,z)\geqslant 0$ are the conduction terms, such that $\sigma_{2}$ contains real physical conductivity as well as additional artificial conductivity related to nonreflective boundary conditions (in this work we use stretched coordinate formulation of the perfectly matched layers, PML, boundary conditions [20]), whereas $\sigma_{1}$ normally contains artificial PML conductivity values only. The functions $J_{H,E}=J_{H,E}(x,y,z,t)$ represent given source terms. If, for the moment, we assume that the homogeneous Dirichlet boundary conditions are supplied with (1), then a standard Yee finite difference discretization on a staggered Cartesian mesh results in a time-continuous space-discretized system

[TABLE]

where the vector functions $\bm{h}(t)$ and $\bm{e}(t)$ contain the mesh values of the unknown fields, $M_{\mu}$, $M_{\varepsilon}$, and $M_{\sigma_{1,2}}$ are diagonal matrices containing the values of $\mu$, $\varepsilon$, and $\sigma_{1,2}$, respectively, $K$ and $K^{T}$ are discrete curl operators and $\bm{j}_{H,E}(t)$ are the mesh values of the source functions $J_{H,E}$.

Note that a semidiscrete system of ordinary differential equations (ODEs), which is very similar to (2), is also obtained when the standard Whitney-Nédélec vector finite elements are employed (see e.g. [26, 7, 31]). In this case $M_{\mu}$, $M_{\varepsilon}$, and $M_{\sigma_{1,2}}$ are the mass matrices. It is convenient to rewrite the system (2) as

[TABLE]

where $M_{1}=M_{\mu}^{-1}M_{\sigma_{1}}$, $K_{1}=M_{\mu}^{-1}K$, $M_{2}=M_{\varepsilon}^{-1}M_{\sigma_{2}}$, $K_{2}^{T}=M_{\varepsilon}^{-1}K^{T}$, and the inverse mass matrices are computed explicitly only if they are diagonal or block diagonal. The latter is the case if discontinuous Galerkin finite elements are used, see e.g. [28]. We denote the size of the ODE system in (3) by $n$, and let $n=n_{1}+n_{2}$, where $n_{1,2}$ are the numbers of degrees of freedom associated with magnetic and electric fields, respectively.

Employment of the nonreflective PML boundary conditions [20, 4] means that auxiliary variables are added to the Maxwell system (1) which, after space discretization, enter the semidiscrete system (3) as well. These additional variables are nonzero only in the so-called PML region (a region just outside the boundary of the domain of interest). Incorporation of the PML boundary conditions into (1),(3) (for a detailed derivation we refer to [8]) leads to the resulting semi-discrete ODE system of an extended size

[TABLE]

where $N=m+n$, with $m$ being the number of space-discretized auxiliary PML variables ($m$ is proportional to the number of mesh points in the PML region), and the matrices $B_{1,2}\in\mathbb{R}^{m\times n}$ couple these variables to the main variables $\bm{h}$ and $\bm{e}$. For representative values of $m$ and $n$ see Table 1 in the numerical experiment section below. The matrices $B_{1,2}$ are defined in more detail below in A.

3 Nested Schur complement solver

Exponential time integration based on rational shift-and-invert Krylov subspace methods [30, 5] as well as implicit time integration [31] of systems (4) involves solution of linear systems

[TABLE]

where $\mathcal{A}$ is defined in (4) and $\gamma>0$ is a given parameter, related (or equal) to the time step size. Matrices having nested saddle point structure111Formally speaking, $I+\gamma\mathcal{A}$ gets a saddle point structure when we switch the sign in the second block row of a linear system $(I+\gamma\mathcal{A})x=b$. as the matrix $\mathcal{A}$ are recently called double saddle point problems [2]. Our starting point in construction of preconditioners for matrices of this type is the observation that an efficient preconditioner should involve a Schur complement (see e.g. [24, 3, 14]). In particular, for matrices

[TABLE]

good Schur complement-based block-diagonal preconditioners are

[TABLE]

For modern Krylov subspace methods such as GMRES, these preconditioners guarantee convergence in three iterations [24]. Applying them for our matrix $I+\gamma\mathcal{A}$ means that linear systems with

[TABLE]

have to be solved efficiently. Comparing these two possible options, we choose option 2 because of the simpler structure of the matrix. Furthermore, assuming for the moment that the systems with the matrix $I+\gamma A$ can be solved efficiently and taking into account that $B_{1}^{T}B_{2}$ is of a low rank, we may expect that $I+\gamma A$ can be a good preconditioner when solving systems with the matrix $I+\gamma A+\gamma^{2}B_{1}^{T}B_{2}$. This expectation is confirmed in practice: number of preconditioned by $I+\gamma A$ iterations to solve systems with $I+\gamma A+\gamma^{2}B_{1}^{T}B_{2}$ remain approximately constant as the spatial discretization mesh gets finer (see Table 2). Moreover, as shown by formula (13) in the appendix below, the matrix $\gamma^{2}B_{1}^{T}B_{2}$ depends on the mesh size in a way similar to the matrix $I+\gamma A$ does: only the $(1,2)$ and $(2,1)$ blocks in this matrix depend on the mesh size as $\sim 1/h$ (assuming $h=h_{x}=h_{y}=h_{z}$).

Now a question arises whether and how the systems with the matrix

[TABLE]

can be solved efficiently. We proceed in a similar way as for the matrix $I+\gamma\mathcal{A}$ and explore the two possible options for a Schur complement-based preconditioner:

[TABLE]

In many applications involving the Maxwell equations (such as, e.g., photonics and gas-and-oil exploration) the permeability $\mu$ is usually constant ($\mu\equiv 1$), whereas $\varepsilon$ is not and can be a strongly varying function. Since the motivation for this work is photonics modeling, we choose for option 1, where the matrix $\gamma^{2}K_{2}^{T}(I+\gamma M_{1})^{-1}K_{1}$ has a simpler structure than the matrix $\gamma^{2}K_{1}(I+\gamma M_{2})^{-1}K_{2}^{T}$ in option 2. It is convenient to rewrite the chosen Schur complement (7) in the form

[TABLE]

which has an advantage that the matrix in brackets is symmetric positive definite. Further inspection of the bracketed matrix reveals that it is similar in structure to a shifted Laplacian, where the shift is given by $M_{\varepsilon}+\gamma M_{\sigma_{2}}$. For this reason a large variety of solvers is available for solving systems with this matrix. These include (i) sparse direct factorization solvers (on coarse meshes); (ii) multigrid solvers (which should not be too difficult in implementation since only one field is involved); (iii) algebraic multigrid methods; (iv) preconditioned conjugate gradients (CG). In this paper we use the CG solver preconditioned by the incomplete Cholesky IC(0) preconditioner [22].

The described nested Schur complement approach can be used either as a Schur complement based preconditioner or as a “direct” solver, computing the inverse action $(I+\gamma\mathcal{A})^{-1}$ with the inner two-level iterative solver for the Schur complement. In Figure 1 we outline the outer part of the introduced nested Schur complement solver for the case the action of $(I+\gamma\mathcal{A})^{-1}$ is computed. The inner part, corresponding to the action of $(I+\gamma A)^{-1}$, can then be computed in a similar fashion, as described above.

4 Other possible preconditioners

A variety of other different preconditioners are available for solving saddle point problems of type (5), see e.g. [3, 14] and recent work [2]. However, a general problem with linear solvers employed in implicit and exponential time integrators is that the additional computational work spent for solving linear systems has to be paid off by an increase in a time step. Assume, for instance, that approximately ten iterations with a basic preconditioner have to be done per time step, such that costs of a preconditioned matrix–vector product (matvec) are approximately equal to the costs of an unpreconditioned matvec (which can be achieved by the Eisenstat’s trick [13]). Then the time step size have to be increased at least by a factor of ten to compensate for the increased costs. Such a time step increase, however, is not always possible due to accuracy restrictions, especially for mildly stiff problems such as the Maxwell equations with PML boundary conditions. This makes the choice of a proper preconditioner difficult and significantly restricts a variety of possible options [31].

Recently an efficient alternating direction implicit (ADI) preconditioner is proposed and analyzed for solving time-dependent Maxwell equations [17] discretized in space by finite differences. In [8] de Cloet, Marissen and Westendorp compared performance of this ADI preconditioner with another preconditioner based on the field splitting (i.e., a splitting into the magnetic and electric fields). Their conclusion is that this field splitting (FS) preconditioner outperforms the ADI preconditioner in terms of the CPU time. Unlike the ADI preconditioner, the FS preconditioner is not restricted to the finite difference approximations on Cartesian meshes. The linear system (5), with either the FS or ADI preconditioner applied from the right, can be written as

[TABLE]

where

[TABLE]

In the FS preconditioner the matrices $\mathcal{A}_{1,2}$ are defined as

[TABLE]

where the matrices $B_{j,H}$, $B_{j,E}$, $j=1,2$, form splittings of the PML blocks $B_{1,2}$,

[TABLE]

defined in B. For complete definition of the ADI preconditioner we refer to [8, Section 5.1] and [17]. In both the ADI and FS preconditioners the factors $I+\gamma\mathcal{A}_{1,2}$ are not triangular matrices (see Figure 2) and sparse LU factorizations have to be carried out to implement the preconditioner actions.

Numerical tests presented in [8] demonstrate that both ADI and FS preconditioners require approximately the same number of iterations to converge. However, the FS preconditioner is faster than ADI in terms of the CPU time. This is caused by another attractive property of the FS preconditioner observed in [8]: sparse LU factorizations of the FS factors $I+\gamma\mathcal{A}_{1,2}$ do not yield any additional fill in. This is the case for the ADI preconditioner, see [8, Section 6.2], where the fill in varies from 25% on coarse meshes to 65% on the mesh $80\times 80\times 48$. Ritz values obtained after 24 iterations of the FOM iterative method [27] with both the ADI and FS preconditioners are plot in Figure 3. There we see that both preconditioners yield preconditioned matrices with effectively real spectra.

5 Numerical experiments

In the test problem a 3D photonic crystal is considered. At the $x$- and $y$- boundaries of the spatial domain $[1,4]\times[1,4]\times[0,3]$ the PML boundary conditions are imposed, whereas homogeneous Dirichlet (perfectly conducting) boundary conditions are posed on the $z$-boundaries. The PML regions extend the total computational domain along the $x$- and $y$-walls to $[0,5]\times[0,5]\times[0,3]$. The crystal consists of $3\times 3\times 3$ spheres of radius $0.4$ centered at points $(x_{i},y_{j},z_{k})=(2.5+i,2.5+j,1.5+k)$, $i,j,k=-1,0,1$. The magnetic permeability $\mu\equiv 1$ in the whole domain, whereas the electric permittivity $\varepsilon$ is set to $8.9$ inside the spheres and to $1$ everywhere else in the domain.

We consider matrices $\mathcal{A}$ resulting from the standard Yee finite difference approximation, see Table 1 for typical mesh sizes used in the tests. The parameter $\gamma$ is chosen as explained in [5] and set to $\gamma=0.012$ in all the tests. The tests are run in Matlab on a Linux PC with 8 CPUs Intel Xeon E5504 2.00GHz.

In Table 2 we illustrate the fact that linear systems with the matrix $I+\gamma A+\gamma^{2}B_{1}^{T}B_{2}$ can be efficiently solved iteratively using $I+\gamma A$ as a preconditioner. The preconditioner actions here are carried with the help of the sparse LU factorization (UMFPACK in Matlab), that is why in this case we can not use a fine mesh. In this test the exact solution vectors are taken to have normally distributed random entries with zero mean and variance one. This is done to make the test difficult so that the solver can not possibly profit from the solution smoothness. The number of iterations listed there are for the BiCGstab(2) iterative solver (the standard built-in solver in Matlab) run to satisfy the stopping criterion tolerance of $10^{-6}$. In the same table we also list the norms of the symmetric and skew-symmetric parts of $\gamma^{2}B_{1}^{T}B_{2}$. The values show that the field of values of this matrix is bounded and confirm the expectation given by relation (13): only the off-diagonal blocks (related to the skew-symmetric part) of this matrix depend on the mesh size and this dependence is linear. In Figure 4 we plot 24 Ritz values of the preconditioned matrix $I+(I+\gamma A)^{-1}\gamma^{2}B_{1}^{T}B_{2}$. As we see, the Ritz values are real and well clustered, which means that the preconditioner is efficient and damps the skew-symmetric part of the system matrix well.

We now test our nested Schur complement solver. The linear systems in the test have the exact solutions which are again a normally distributed random vector whose entries have zero mean and variance one. This is done to prevent the solver from profiting from the solution smoothness. The solver is applied in its direct form, as outlined in Figure 1, with the outer solver GMRES(10) and inner solver ICCG(0). The stopping criterion tolerance in both solvers is set to $10^{-10}$. We see that the number of the outer iterations remains constant, independently of the mesh size, as expected. The number of inner iterations grows because the CG solver is used with the simple IC(0) preconditioner. We note that the number of inner iterations changes from one outer iteration to another and the inner iteration count reported in the table is the maximum number of inner iterations (required in all cases at the last outer iteration).

Finally, in Table 4 we present results of comparison of the nested Schur complement solver with the FS preconditioner. The comparison is done on linear systems arising in the time integration carried out by an exponential integrator based on a rational shift-and-invert exponential Krylov subspace method. Therefore the stopping criterion tolerance in the tests vary and is reported in the table. We see that the nested Schur complement solver outperforms the FS preconditioner on fine meshes. This is expected because the FS preconditioner does not converge mesh-independently.

6 Conclusions and an outlook to further research

A nested Schur complement solver is proposed for iterative linear system solution within exponential and implicit time integration of the Maxwell equations. The solver exhibits a mesh-independent convergence and outperforms other preconditioners, such ADI (alternative direction implicit) and FS (field splitting) preconditioners.

Different aspects in the proposed concept require further investigation and possible improvement. In future we plan to test the nested Schur complement solver in combination with a more robust (and mesh-independent) inner iterative solver. Another interesting research question is which form of the solver is most efficient: its direct form, as tested in this paper, or iterative, as a three-level iterative solver.

Appendix A PML matrices $B_{1,2}$

For simplicity of the presentation we assume that the Yee finite difference discretization is used on a mesh of $n_{x}\times n_{y}\times n_{z}$ cells. Then it is possible to construct this discretization in such a way that all the three components of each of the two fields have the same number of degrees of freedom $(n_{x}+1)(n_{y}+1)(n_{z}+1)$, This can be done by augmenting the “shorter” components with auxiliary “void” degrees of freedom. Then $n_{1}=n_{2}=n/2$ and the matrices $B_{1,2}$ are defined with the help of the following matrices $\widehat{B}_{1,2}\in\mathbb{R}^{2n\times n}$:

[TABLE]

where diagonal matrices $\Sigma,\Sigma_{*}\in\mathbb{R}^{n_{1}\times n_{1}}$ read

[TABLE]

Here $\bm{\sigma}_{x,y,z}$ are the PML artificial conductivity values for the corresponding direction and $\cdot$ denotes elementwise multiplication. Since $\bm{\sigma}_{x,y,z}$ are only nonzero inside the PML region, there are a lot of redundant zero rows in $\widehat{B}_{1,2}$. By omitting these zero rows we obtain the matrices $B_{1,2}$. Note also that

[TABLE]

which means that only the $(1,2)$ and $(2,1)$ blocks in this matrix depend on the mesh (as $K_{1,2}$ are discrete curl operators) as $\sim 1/h$ (assuming $h=h_{x}=h_{y}=h_{z}$). This is important in construction of preconditioners for linear system solution with the matrix $I+\gamma A+\gamma^{2}B_{1}^{T}B_{2}$.

Appendix B Field splitting of the matrices $B_{1,2}$

To define the FS preconditioner we split the matrices $\widehat{B}_{1,2}$ into the components corresponding to the two fields $H$ and $E$ as follows:

[TABLE]

where

[TABLE]

and $\widehat{B}_{1,E}:=\widehat{B}_{1}-\widehat{B}_{1,H}$, $\widehat{B}_{2,E}:=\widehat{B}_{2}-\widehat{B}_{2,H}$. The matrices ${B}_{1,H}$, ${B}_{2,H}$, ${B}_{1,E}$, ${B}_{2,E}$ are then obtained from the matrices $\widehat{B}_{1,H}$, $\widehat{B}_{2,H}$, $\widehat{B}_{1,E}$, $\widehat{B}_{2,E}$, respectively, by omitting zero rows in the latter four.