Multitrace formulations and Domain Decomposition Methods for the solution of Helmholtz transmission problems for bounded composite scatterers

TL;DR

This paper develops and analyzes multitrace formulations and domain decomposition methods for efficiently solving Helmholtz transmission problems in composite scatterers, demonstrating improved iterative solver performance.

Contribution

It introduces a novel DDM approach with generalized Robin boundary conditions using Fourier multiplier approximations, enhancing solver efficiency for Helmholtz problems.

Findings

Hierarchical elimination improves Krylov solver efficiency

Generalized Robin conditions reduce iteration counts

Method scales mildly with frequency and subdomain number

Abstract

We present Nystr\"om discretizations of multitrace formulations and non-overlapping Domain Decomposition Methods (DDM) for the solution of Helmholtz transmission problems for bounded composite scatterers with piecewise constant material properties. We investigate the performance of DDM with both classical Robin and generalized Robin boundary conditions. The generalized Robin boundary conditions incorporate square root Fourier multiplier approximations of Dirichlet to Neumann operators. While the classical version of DDM is not particularly well suited for Krylov subspace iterative solvers, we show that the associated DDM linear system can be efficiently solved by hierarchical elimination via Schur complements of the Robin data. We show through numerical examples that the latter version of DDM gives rise to small numbers of Krylov subspace iterations that depend mildly on the frequency and number of subdomains.

Figures

Tables

Table 1. Table 1: Performance of the various formulations considered in this text in the three subdomain case with ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 64 subscript 𝜀 1 64 \varepsilon_{1}=64 , and ε 2 = 256 subscript 𝜀 2 256 \varepsilon_{2}=256 . The numbers of unknowns required by the DDM and g-MTF are 384 384 384 , 768 768 768 , and 1 , 536 1 536 1,536 respectively; the MTF uses twice as many unknowns in each case. The number of terms used in the Neumannn series ( 4.13 ) were N t ​ r ​ u ​ n ​ c = 40 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 40 N_{trunc}=40 , N t ​ r ​ u ​ n ​ c = 80 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 80 N_{trunc}=80 , and N t ​ r ​ u ​ n ​ c = 80 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 80 N_{trunc}=80 in the cases ω = 1 𝜔 1 \omega=1 , ω = 2 𝜔 2 \omega=2 , and ω = 4 𝜔 4 \omega=4 respectively.

$\omega$	DDM (4.6)				DDM Schur (4.10)		DDM Schur (4.13)		MTF				g-MTF (3.2)
	It	It DtNR	It GSqR	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It (3.8)	It (6.1)	It (6.2)	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$
1	137	16	25	4.8 $\times$ ${10}^{-3}$	15	4.9 $\times$ ${10}^{-3}$	15	3.2 $\times$ ${10}^{-2}$	431	242	160	5.3 $\times$ ${10}^{-3}$	159	3.9$\times {10}^{-3}$
2	182	19	40	1.8 $\times$ ${10}^{-3}$	18	2.6 $\times$ ${10}^{-3}$	19	4.8 $\times$ ${10}^{-2}$	750	433	293	2.4 $\times$ ${10}^{-3}$	270	3.2$\times {10}^{-3}$
4	282	19	68	2.9 $\times$ ${10}^{-3}$	20	1.3 $\times$ ${10}^{-3}$	21	1.0 $\times$ ${10}^{-2}$	1,335	732	508	3.4 $\times$ ${10}^{-3}$	430	2.8$\times {10}^{-3}$

Table 2. Table 2: Performance of the various formulations considered in this text in the five subdomain case with ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 64 subscript 𝜀 1 64 \varepsilon_{1}=64 , ε 2 = 256 subscript 𝜀 2 256 \varepsilon_{2}=256 , ε 3 = 1024 subscript 𝜀 3 1024 \varepsilon_{3}=1024 , and ε 4 = 4096 subscript 𝜀 4 4096 \varepsilon_{4}=4096 . The numbers of unknowns required by the DDM and g-MTF are 4 , 608 4 608 4,608 and 9 , 216 9 216 9,216 respectively; the MTF uses twice as many unknowns in each case. The number of terms used in the Neumannn series ( 4.13 ) were N t ​ r ​ u ​ n ​ c = 160 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 160 N_{trunc}=160 in both cases ω = 1 𝜔 1 \omega=1 , and ω = 2 𝜔 2 \omega=2 respectively.

$\omega$	DDM (4.6)				DDM Schur (4.10)		DDM Schur (4.13)		MTF				g-MTF (3.2)
	It	It DtNR	It GSqR	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It (3.8)	It (6.1)	It (6.2)	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$
1	977	39	70	1.6 $\times$ ${10}^{-3}$	12	1.9 $\times$ ${10}^{-3}$	14	3.5 $\times$ ${10}^{-2}$	2,419	1,296	1,107	4.5 $\times$ ${10}^{-3}$	942	1.1$\times$ ${10}^{-3}$
2	1,788	72	140	2.9 $\times$ ${10}^{-3}$	20	3.3 $\times$ ${10}^{-3}$	23	2.8 $\times$ ${10}^{-2}$	4,534	2,442	2,154	5.4 $\times$ ${10}^{-3}$	1,782	2.1$\times$ ${10}^{-3}$

Table 3. Table 3: Performance of the various formulations considered in this text in the three subdomain case with ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 4 subscript 𝜀 1 4 \varepsilon_{1}=4 , and ε 2 = 16 subscript 𝜀 2 16 \varepsilon_{2}=16 . The numbers of unknowns required by the DDM and g-MTF are 384 384 384 , 768 768 768 , 1 , 536 1 536 1,536 , and 3 , 072 3 072 3,072 respectively; the MTF uses twice as many unknowns in each case. The number of terms used in the Neumannn series ( 4.13 ) were N t ​ r ​ u ​ n ​ c = 20 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 20 N_{trunc}=20 , N t ​ r ​ u ​ n ​ c = 20 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 20 N_{trunc}=20 , N t ​ r ​ u ​ n ​ c = 30 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 30 N_{trunc}=30 , and N t ​ r ​ u ​ n ​ c = 80 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 80 N_{trunc}=80 in the cases ω = 4 𝜔 4 \omega=4 , ω = 8 𝜔 8 \omega=8 , ω = 16 𝜔 16 \omega=16 , and ω = 32 𝜔 32 \omega=32 respectively.

$\omega$	DDM (4.6)				DDM Schur (4.10)		DDM Schur (4.13)		MTF				g-MTF (3.2)
	It	It DtNR	It GSqR	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It (3.8)	It (6.1)	It (6.2)	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$
4	88	25	27	3.6 $\times$ ${10}^{-3}$	25	3.5 $\times$ ${10}^{-3}$	26	8.5 $\times$ ${10}^{-3}$	212	150	109	4.3 $\times$ ${10}^{-3}$	108	1.8$\times {10}^{-3}$
8	142	40	41	4.4 $\times$ ${10}^{-3}$	40	4.3 $\times$ ${10}^{-3}$	42	4.2 $\times$ ${10}^{-2}$	342	235	176	5.8 $\times$ ${10}^{-3}$	174	4.1$\times {10}^{-3}$
16	176	49	57	3.1 $\times$ ${10}^{-3}$	45	3.0 $\times$ ${10}^{-3}$	48	4.9 $\times$ ${10}^{-2}$	483	361	267	6.3 $\times$ ${10}^{-3}$	280	5.8$\times {10}^{-3}$
32	274	64	91	1.9 $\times$ ${10}^{-3}$	59	1.3 $\times$ ${10}^{-3}$	64	4.5 $\times$ ${10}^{-2}$	808	565	405	3.8 $\times$ ${10}^{-3}$	475	7.8$\times {10}^{-3}$

Table 4. Table 4: Performance of the various formulations considered in this text in the five subdomain case with ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 4 subscript 𝜀 1 4 \varepsilon_{1}=4 , ε 2 = 16 subscript 𝜀 2 16 \varepsilon_{2}=16 , ε 3 = 64 subscript 𝜀 3 64 \varepsilon_{3}=64 , and ε 4 = 256 subscript 𝜀 4 256 \varepsilon_{4}=256 . The numbers of unknowns required by the DDM and g-MTF are 1 , 152 1 152 1,152 , 2 , 304 2 304 2,304 and 4 , 608 4 608 4,608 respectively; the MTF uses twice as many unknowns in each case. The number of terms used in the Neumannn series ( 4.13 ) were N t ​ r ​ u ​ n ​ c = 80 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 80 N_{trunc}=80 , N t ​ r ​ u ​ n ​ c = 80 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 80 N_{trunc}=80 , and N t ​ r ​ u ​ n ​ c = 160 subscript 𝑁 𝑡 𝑟 𝑢 𝑛 𝑐 160 N_{trunc}=160 in the cases ω = 4 𝜔 4 \omega=4 , ω = 8 𝜔 8 \omega=8 , and ω = 16 𝜔 16 \omega=16 respectively.

$\omega$	DDM (4.6)				DDM Schur (4.10)		DDM Schur (4.13)		MTF				g-MTF (3.2)
	It	It DtNR	It GSqR	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$	It (3.8)	It (6.1)	It (6.2)	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$
4	427	52	67	2.3 $\times$ ${10}^{-3}$	23	2.8 $\times$ ${10}^{-3}$	30	3.6 $\times$ ${10}^{-2}$	1,152	673	599	6.0 $\times$ ${10}^{-3}$	859	5.4$\times {10}^{-3}$
8	613	80	109	1.1 $\times$ ${10}^{-3}$	32	1.3 $\times$ ${10}^{-3}$	40	4.5 $\times$ ${10}^{-2}$	2,058	1,227	1,111	6.7 $\times$ ${10}^{-3}$	1,426	4.9$\times {10}^{-3}$
16	784	106	178	9.3 $\times$ ${10}^{-4}$	41	2.1 $\times$ ${10}^{-3}$	51	3.1 $\times$ ${10}^{-2}$	3,371	2,163	1,997	5.3 $\times$ ${10}^{-3}$	2,471	5.4$\times {10}^{-3}$

Table 5. Table 5: Performance of the various formulations considered in this text in the three subdomain case with ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 4 subscript 𝜀 1 4 \varepsilon_{1}=4 , ε 2 = 16 subscript 𝜀 2 16 \varepsilon_{2}=16 and ω = 2 𝜔 2 \omega=2 . The numbers of unknowns required by the DDM and g-MTF are 192 192 192 and 384 384 384 respectively; the MTF uses twice as many unknowns in each case. We present results for the DDM formulations in the case when the subdomains Ω 1 subscript Ω 1 \Omega_{1} and Ω 2 subscript Ω 2 \Omega_{2} are each subdivided into two subdomains, and we refer to these cases by 4 DDM.

DDM (4.6)		DDM Schur (4.10)	4 DDM (4.6)		4 DDM Schur (4.10)	MTF				g-MTF (3.2)
It	${\epsilon }_{\infty }$	It	It	${\epsilon }_{\infty }$	It	It (3.8)	It (6.1)	It (6.2)	${\epsilon }_{\infty }$	It	${\epsilon }_{\infty }$
77	4.1 $\times$ ${10}^{-3}$	22	107	4.5 $\times$ ${10}^{-3}$	23	142	102	77	3.6 $\times$ ${10}^{-3}$	67	1.4$\times {10}^{-3}$
73	3.4 $\times$ ${10}^{-4}$	22	118	3.5 $\times$ ${10}^{-4}$	22	145	104	78	2.8 $\times$ ${10}^{-4}$	69	1.7$\times {10}^{-4}$

Table 6. Table 6: Left panel: numbers of GMRES iterations needed to reach a residual of 10 − 4 superscript 10 4 10^{-4} for the solution of equation ( 5 ) in each subdomain and respectively for the solution of the global DDM system ( 4.3 ) in the three subdomain case with ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 64 subscript 𝜀 1 64 \varepsilon_{1}=64 , and ε 2 = 256 subscript 𝜀 2 256 \varepsilon_{2}=256 (top three rows) and ε 0 = 1 subscript 𝜀 0 1 \varepsilon_{0}=1 , ε 1 = 4 subscript 𝜀 1 4 \varepsilon_{1}=4 , and ε 2 = 16 subscript 𝜀 2 16 \varepsilon_{2}=16 (bottom four rows). Right panel: same numbers of GMRES iterations when the subdomains Ω 1 subscript Ω 1 \Omega_{1} and Ω 2 subscript Ω 2 \Omega_{2} are each subdivided into two half-sized non-overlapping subdomains Ω 1 = Ω 11 ∪ Ω 12 subscript Ω 1 subscript Ω 11 subscript Ω 12 \Omega_{1}=\Omega_{11}\cup\Omega_{12} and Ω 2 = Ω 21 ∪ Ω 22 subscript Ω 2 subscript Ω 21 subscript Ω 22 \Omega_{2}=\Omega_{21}\cup\Omega_{22} respectively.

$\omega$	DDM (4.6)						4 DDM (4.6)
	It ${\Omega }_0$	It ${\Omega }_1$	It ${\Omega }_2$	It	It DtNR	It GSqR	It ${\Omega }_0$	It ${\Omega }_{11}{\Omega }_{12}$	It ${\Omega }_{21}{\Omega }_{22}$	It	It DtNR	It GSqR
1	8	30	48	137	16	25	8	27	39	225	42	44
2	9	47	79	182	19	40	9	39	64	361	61	73
4	11	70	121	282	19	68	10	59	109	597	104	129
4	10	21	36	88	31	30	10	19	32	153	46	37
8	13	25	46	142	40	54	12	21	40	231	64	63
16	16	32	60	176	49	77	15	31	56	257	78	78
32	19	39	73	274	64	106	18	38	70	367	101	130

Equations

$$\Delta u_j+k_j^2{u}_j$$

$$u_j+{\delta }_j^0{u}^{inc}$$

$${\epsilon }_j^{-1}({\partial }_{n_j}{u}_j+{\delta }_j^0{\partial }_{n_j}{u}^{inc})$$

$$\underset{r\to \infty }{\lim }{r}^{1/2}(\partial u_0/\partial r-i{k}_0{u}_0)$$

$$[S{L}_{\Gamma ,k}(\phi )](\mathbf{z}):={\int }_{\Gamma }{G}_k(\mathbf{z}-\mathbf{y})\phi (\mathbf{y})ds(\mathbf{y}),\mathbf{z}\in {\mathbb{R}}^2\setminus \Gamma ,$$

$$[D{L}_{\Gamma ,k}(\phi )](\mathbf{z}):={\int }_{\Gamma }\frac{\partial G_k(\mathbf{z}-\mathbf{y})}{\partial \mathbf{n}(\mathbf{y})}\phi (\mathbf{y})ds(\mathbf{y}),\mathbf{z}\in {\mathbb{R}}^2\setminus \Gamma ,$$

$${\gamma }_{\Gamma }^{D,1}S{L}_{\Gamma ,k}(\phi )$$

$${\gamma }_{\Gamma }^{N,j}S{L}_{\Gamma ,k}(\phi )$$

$$u_j=-D{L}_{\partial {\Omega }_j,j}{u}_j+S{L}_{\partial {\Omega }_j,j}{\partial }_{n_j}{u}_j\text{mbox}in{\Omega }_j.$$

$$\frac{1}{2}{u}_0+K_{\partial {\Omega }_0,0}{u}_0-{\epsilon }_0{S}_{\partial {\Omega }_0,0}({\epsilon }_0^{-1}{\partial }_{n_0}{u}_0)=0\mathrm{on}\partial {\Omega }_0.$$

$$u_0=\{\begin{array}{ll}{u}_1-u^{inc},& \mathrm{on}{\Gamma }_{01},\\ u_2-u^{inc},& \mathrm{on}{\Gamma }_{02},\end{array}$$

$$u_0=X_{01}{u}_1+X_{02}{u}_2-u^{inc},\mathrm{on}\partial {\Omega }_0,$$

$$X_{01}:=E_{01}^{(0)}{\Pi }_{01}{X}_{02}:=E_{02}^{(0)}{\Pi }_{02}$$

$$E_{01}^{(0)}{\psi }_{01}:=\{\begin{array}{ll}{\psi }_{01},& \mathrm{on}{\Gamma }_{01}\\ 0,& \mathrm{on}{\Gamma }_{02}\end{array}{E}_{02}^{(0)}{\psi }_{02}:=\{\begin{array}{ll}0,& \mathrm{on}{\Gamma }_{01}\\ {\psi }_{02},& \mathrm{on}{\Gamma }_{02}\end{array}$$

$$K_{\partial {\Omega }_0,0}{u}_0-{\epsilon }_0{S}_{\partial {\Omega }_0,0}({\epsilon }_0^{-1}{\partial }_{n_0}{u}_0)+\frac{1}{2}{X}_{01}{u}_1+\frac{1}{2}{X}_{02}{u}_2=\frac{1}{2}{u}^{inc}\mathrm{on}\partial {\Omega }_0.$$

$$\frac{1}{2}{\epsilon }_0^{-1}{\partial }_{n_0}{u}_0-K_{\partial {\Omega }_0,0}^{\top }({\epsilon }_0^{-1}{\partial }_{n_0}{u}_0)+{\epsilon }_0^{-1}{N}_{\partial {\Omega }_0,0}{u}_0=0\mathrm{on}\partial {\Omega }_0.$$

$${\epsilon }_0^{-1}{\partial }_{n_0}{u}_0:=\{\begin{array}{ll}-{\epsilon }_1^{-1}{\partial }_{n_1}{u}_1-{\epsilon }_0^{-1}{\partial }_{n_0}{u}^{inc},& \mathrm{on}{\Gamma }_{01},\\ -{\epsilon }_2^{-1}{\partial }_{n_2}{u}_2-{\epsilon }_0^{-1}{\partial }_{n_0}{u}^{inc},& \mathrm{on}{\Gamma }_{02},\end{array}$$

$${\epsilon }_0^{-1}{\partial }_{n_0}{u}_0=-X_{01}({\epsilon }_1^{-1}{\partial }_{n_1}{u}_1)-X_{02}({\epsilon }_2^{-1}{\partial }_{n_2}{u}_2)-{\epsilon }_0^{-1}{\partial }_{n_0}{u}^{inc},\mathrm{on}\partial {\Omega }_0,$$

$${\epsilon }_0^{-1}{N}_{\partial {\Omega }_0,0}{u}_0-K_{\partial {\Omega }_0,0}^{\top }({\epsilon }_0^{-1}{\partial }_{n_0}{u}_0)-\frac{1}{2}{X}_{01}({\epsilon }_1^{-1}{\partial }_{n_1}{u}_1)-\frac{1}{2}{X}_{02}({\epsilon }_2^{-1}{\partial }_{n_2}{u}_2)=\frac{1}{2}{\epsilon }_0^{-1}{\partial }_{n_0}{u}^{inc}\mathrm{on}\partial {\Omega }_0.$$

$${\mathcal{C}}_j:=\left[\begin{array}{cc}{K}_{\partial {\Omega }_j,j}& -{\epsilon }_j{S}_{\partial {\Omega }_j,j}\\ {\epsilon }_j^{-1}{N}_{\partial {\Omega }_j,j}& -K_{\partial {\Omega }_j,j}^{\top }\end{array}\right],j=0,1,2,{\mathcal{C}}_j^2=\frac{1}{4}\left[\begin{array}{cc}I& 0\\ 0& I\end{array}\right],$$

$${\mathbb{X}}_{01}=\left[\begin{array}{cc}{X}_{01}& 0\\ 0& -X_{01}\end{array}\right]{\mathbb{X}}_{02}=\left[\begin{array}{cc}{X}_{02}& 0\\ 0& -X_{02}\end{array}\right].$$

$${\mathbf{u}}_j:=\left[\begin{array}{c}{u}_j\\ {\epsilon }_j^{-1}{\partial }_{n_j}{u}_j\end{array}\right],j=0,1,2,$$

$${\mathbf{u}}^{inc}:=\left[\begin{array}{c}{u}^{inc}\\ {\epsilon }_0^{-1}{\partial }_{n_0}{u}^{inc}\end{array}\right],$$

$${\mathcal{C}}_0{\mathbf{u}}_0+\frac{1}{2}{\mathbb{X}}_{01}{\mathbf{u}}_1+\frac{1}{2}{\mathbb{X}}_{02}{\mathbf{u}}_2=\frac{1}{2}{\mathbf{u}}^{inc}\mathrm{on}\partial {\Omega }_0.$$

$$\left[\begin{array}{ccc}{\mathcal{C}}_0& {\mathbb{X}}_{01}& {\mathbb{X}}_{02}\\ {\mathbb{X}}_{10}& {\mathcal{C}}_1& {\mathbb{X}}_{12}\\ {\mathbb{X}}_{20}& {\mathbb{X}}_{21}& {\mathcal{C}}_2\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_0\\ {\mathbf{u}}_1\\ {\mathbf{u}}_2\end{array}\right]=\frac{1}{2}\left[\begin{array}{c}{\mathbf{u}}^{inc}\\ -{\mathbb{X}}_{10}{\mathbf{u}}^{inc}\\ -{\mathbb{X}}_{20}{\mathbf{u}}^{inc}\end{array}\right],$$

$$u_j(\mathbf{x}):=S{L}_{\Gamma ,j}v-{\epsilon }_{j}D{L}_{\Gamma ,j}p,\mathbf{x}\in {\Omega }_j$$

$$\frac{{\epsilon }_{\ell }+{\epsilon }_j}{2}p+(-1{)}^{\ell }[S_{\ell }-S_j]v+(-1{)}^{\ell }[{\epsilon }_{\ell }{K}_{\ell }-{\epsilon }_j{K}_j]p$$

$$\frac{{\epsilon }_{\ell }^{-1}+{\epsilon }_j^{-1}}{2}v+(-1{)}^{\ell }[{\epsilon }_j^{-1}{S}_j-{\epsilon }_{\ell }^{-1}{S}_{\ell }]v+(-1{)}^{\ell }[N_j-N_{\ell }]p$$

$$\Delta u_j+k_j^2{u}_j$$

$${\epsilon }_j^{-1}({\partial }_{n_j}{u}_j+{\delta }_j^0{\partial }_{n_j}{u}^{inc})+i\eta (u_j+{\delta }_j^0{u}^{inc})$$

$${\mathcal{S}}^j({\psi }_j):=({\epsilon }_j^{-1}{\partial }_{n_j}{u}_j-i\eta u_j){|}_{\partial {\Omega }_j},j=0,1,2,$$

$$\Delta u_j+k_j^2{u}_j$$

$${\epsilon }_j^{-1}{\partial }_{n_j}{u}_j+i\eta u_j$$

$$f_j:=({\epsilon }_j^{-1}{\partial }_{n_j}{u}_j+i\eta u_j){|}_{\partial {\Omega }_j},j=0,1,2$$

Full text

Keywords: multiple junctions, multitrace formulations, domain decomposition methods.

AMS subject classifications: 65N38, 35J05, 65T40,65F08

1 Introduction

The phenomenon of electromagnetic wave scattering by bounded penetrable objects composed of several subdomains with different but constant electric permittivities is relevant for numerous applications in antenna design, diffraction gratings, and photovoltaic cells, to name but a few. It is typical in all these application areas that multiple media meet at a single point, a scenario that is referred to as triple or multiple junctions. Standard numerical methods perform poorly when dealing with wave scattering by composite objects with piecewise constant material properties portraying large frequency ranges, as they need to resolve wave interactions with high-contrast sharp interfaces. Volumetric discretizations of these problems result in very large linear systems of equations that are ill-conditioned in the high-frequency regime and whose solution by iterative solvers require inordinate numbers of iterations. Several preconditioning strategies have been proposed to mitigate the aforementioned issues, the most successful arguably being those that rely on the shifted Laplacean [3, 14] or the sweeping preconditioner introduced in [13].

Domain Decomposition Methods (DDM) are natural candidates for the solution of scattering problems involving composite scatterers. DDM are divide and conquer strategies whereby the computational domain is divided into smaller subdomains for which solutions are matched via transmission conditions on subdomain interfaces. The convergence of DDM for time-harmonic wave scattering applications depends a great deal on the choice of the transmission conditions that allow the exchange of information between adjacent subdomains. These interface transmission conditions should ideally allow information to flow out of a subdomain with as little as possible information being reflected back into the subdomain. Thus, the interface transmission conditions fall into the category of Absorbing Boundary Conditions (ABC). From this perspective, the ideal choice of transmission conditions on an interface between two subdomains is such that the impedance/transmission operator is the restriction to the common interface of the Dirichlet-to-Neumann (DtN) operator corresponding to the adjacent subdomain. Traditionally, the interface transmission conditions were chosen as the classical first order ABC outgoing Robin/impedance boundary conditions [11, 16]. The convergence of DDM with classical Robin interface boundary conditions is slow and is adversely affected by the number of subdomains. Fortunately, the convergence of DDM can be considerably improved through incorporation of ABC that constitute higher order approximations of DtN operators in the form of second order approximations with optimized tangential derivative coefficients [15], square root approximations [4], or other types of non-local transmission conditions [16, 30]. Alternatively, so-called Perfectly Matched Layers can be used at subdomain interfaces [32]. Although the use of more sophisticated ABC, as those recounted above, accelerates a great deal the convergence of DDM, the number of iterations required for convergence still grows –albeit not drastically– with the frequency and number of subdomains. This is not entirely surprising since higher order ABC only involve local exchange of information between adjacent subdomains, and affect to a lesser degree the global exchange of information between distant subdomains. Recent efforts have been devoted to construct “double sweep”-type preconditioners that address the latter issue [34, 35]. The resulting preconditioned DDM scale favorably with frequency and number of subdomains, but appear to be somewhat less effective for wave propagation problems in composite media that exhibit sharp high-contrast interfaces.

Boundary integral equation based solvers for scattering by composite objects with piecewise constant material properties require significantly fewer unknowns than volumetric solvers as only the interfaces of material discontinuity need be discretized. The formulation of these problems in terms of robust boundary integral equations has recently received significant interest in the community, the main achievement being the introduction of Multi-Trace Formulations (MTF) [19, 8]. The derivation of one of the MTF –dubbed local– consists of the following steps: (1) use of Green’s identities in each subdomain, whose boundary is a union of interfaces of material discontinuity, to represent the fields in that subdomain via layer potential; (2) application of Dirichlet and Neumann traces associated to that subdomain to the Green’s identities; followed by (3) enforcement of the continuity conditions across interfaces to replace the identity terms in the previous steps by Dirichlet and Neumann traces of solutions in adjacent subdomains. This procedure leads to a boundary integral equation of the first kind whose unknowns are both interior and exterior Dirichlet and Neumann traces of fields on each interface and which involves, the standard Boundary Integral Operators (BIO) on each subdomain corresponding to the wavenumber associated with that subdomain. Galerkin discretizations of the local MTF using $h$-refinement order bases for relative low frequencies preconditioned by diagonal or Calderón preconditioners yield small Krylov subspace iteration numbers. However, as contrast ratios increase, the Calderón preconditioners are less effective when used for low order Galerkin discretizations. A remedy is to resort to spectral or high-order Galerkin discretizations [21, 22]. The MTF of the second kind—or global version—can be derived if the fields are sought in terms of suitable linear combinations of layer potentials defined on the union of all interfaces–typically referred to as the skeleton–of material discontinuity [8, 18].

In this paper, we illustrate the performance of Nyström solvers based multitrace formulations and DDM solvers in the case of high-frequency scattering problems from composite high-contrast scatterers. A major advantage of MTF is the ease with which they can be incorporated into existing boundary integral equation solvers. We present in this work a straightforward extension of the Helmholtz transmission Nyström solvers introduced in [12] to MTF. We investigate two types of simple preconditioners for the MTF of the first kind: the Calderón preconditioner proposed in [19] and a new preconditioner that combines a Schur complement approach with a Calderón preconditioner. Both preconditioners are shown to be effective for high-frequency high-contrast scattering problems from composite scatterers not depending on the number of subdomains. However, the numbers of iterations required by Nyström discretizations of MTF grows considerably with the frequency and/or the contrast between subdomains, even after resorting to preconditioning. We show that in the aforementioned frequency regime DDM based on boundary integral equations can be advantageous alternatives to MTF. We investigate both DDM based on the exchange of classical Robin data between subdomain as well as DDM that incorporate general Robin/impedance boundary operators that are square root approximations of DtN maps [20]. In both cases, we solve the subdomain Helmholtz equations with Robin or generalized Robin boundary conditions using well-conditioned boundary integral formulations solved by Nyström discretizations. Provided the size of the subdomains is small enough–in terms of wavelengths across–, the latter problems can be solved by direct linear algebra methods.

The numerical results presented in this paper corroborate the by now well-known fact that DDM based on classical Robin boundary conditions do not perform well as either the frequency or/and the number of subdomains are increased. We present two possible approaches to overcome these shortcomings. The first approach consist of hierarchical elimination from the DDM linear system via Schur complements of the Robin data corresponding to interior subdomain interfaces. This idea was presented recently as a hierarchical merging of Robin-to-Robin (RtR) maps [17] in the context of scattering from variable media with smooth index of refraction, and was shown to be equivalent to hierarchical elimination of interior Robin data from the DDM linear system corresponding to multiple scattering [28]. The advantage of the Schur complement DDM is that at each stage inverses of relatively small matrices need be computed, and these can be performed in a hierarchical fashion that optimizes the computational cost without compromising the inherent parallelism of DDM. The second approach we pursue is to incorporate square root approximations–using complex wavenumbers–of the DtN operators as generalized Robin operators in DDM. The operators were used successfully for regularizing boundary integral equations at high frequencies [2, 1, 6, 5, 12, 7] and also incorporated successfully in DDM [4, 34]. Since for composite scatterers, adjacent domain decomposition subdomains can have different wavenumbers, we blend the approximations of DtN operators corresponding to different subdomains–with different wavenumbers as well–via smooth cutoff functions.

Important advantages of using Nyström discretizations are (a) the square root DtN approximations which are naturally defined in the Fourier domain can be implemented with ease as Fourier multipliers using FFTs and (b) the regularization strategy presented in [33] is applicable to general Robin boundary conditions that involve the square root operators. In contrast, boundary element or other type of discretizations (e.g. finite difference or finite element) require either local approximations or rational Padé approximations of the square root operators, and these approximations additional require tangential derivatives [4, 34]. In addition, neither the presence of cross points at multiple junctions nor the incorporation of square root DtN approximations as generalized Robin boundary conditions affect the accuracy of the DDM solvers. We note that the accuracy of finite difference/finite element DDM based on generalized Robin boundary conditions has not been reported in the literature [4, 34]. The Nyström DDM with generalized Robin boundary conditions requires much fewer iterations than the classical Robin boundary conditions DDM while the computational cost of the local subdomain solves for the former method is virtually the same as that for the latter.

The structure of this paper is as follows. Section 2 describes the problems of scattering considered in this paper. Section 2 reviews both first- and second-kind MTF approaches. The proposed DDM approach is then introduced and analyzed in Section 3. Section 5, on the other hand, discusses high-order Nytström discretizations for the various integral equation formulations considered, while a variety of numerical results are shown in Section 6. Finally, the conclusions of this work are presented in Section 7.

2 Electromagnetic transmission problems in composite domains

We consider the problem of two dimensional electromagnetic scattering by structures that feature multiple junctions, i.e. points where more than three interfaces of material discontinuity meet (e.g. the structure displayed in Figure 1). For the sake of presentation simplicity, we focus our treatment of transmission problems with multiple junction domains on the two subdomain case depicted in Figure 1. Specifically, we seek to solve the Helmholtz transmission problem that consists of finding the fields $u_{0}$, $u_{1}$, and $u_{2}$ solutions of the system of equations:

[TABLE]

where $\delta_{j}^{0}$ and $\delta_{\ell}^{0}$ stand for Kronecker operators, that is $\delta_{j}^{0}$ is the identity operator if $j=0$ and the null operator otherwise. Here, the wavenumber $k_{j}$ in the subdomain $\Omega_{j}$ is given by $k_{j}=\omega\sqrt{\varepsilon_{j}}$ in terms of the angular frequency $\omega>0$ and the electric permittivity $\varepsilon_{j}$. Throughout this paper we assume that all the permittivities $\varepsilon_{j}$ are positive real numbers, extensions to more general cases being straightforward. The unit normal to the boundary $\partial\Omega_{j}$ is here denoted by $n_{j}$ and is assumed to point to the exterior of the subdomain $\Omega_{j}$. The incident field $u^{inc}$, on the other hand, is assumed to satisfy the Helmholtz equation with wavenumber $k_{0}$ in the unbounded domain $\Omega_{0}$. Finally, we point out that the well-posedness of the transmission problem (2.1)is discussed in [19] and references therein.

In what follows, we review three main formulations of the transmission problem (2.1). Two of them rely on boundary integral equations, and the third is a domain decomposition method.

3 Multitrace formulations

In this section we review the derivation of the MTF [19] of the transmission problem (2.1). To this end, we make use of the four BIO associated with the Calderón calculus for a Lipschitz domain. Let $D$ be a bounded domain in $\mathbb{R}^{2}$ whose boundary $\Gamma$ is a closed Lipschitz curve. Given a wavenumber $k$, and a density $\varphi$ defined on $\Gamma$, we recall the definitions of the single layer potential [26]

[TABLE]

and the double layer potential

[TABLE]

where $G_{k}(\mathbf{x})=\frac{i}{4}H_{0}^{(1)}(k|\mathbf{x}|)$ represents the two-dimensional Green’s function of the Helmholtz equation with wavenumber $k$, and $\mathbf{n}$ represents the unit normal pointing outside the domain $D$. Applying Dirichlet and Neumann exterior and interior traces on $\Gamma$–denoted by $\gamma_{\Gamma}^{D,1}$ and $\gamma_{\Gamma}^{D,2}$ and respectively $\gamma_{\Gamma}^{N,1}$ and $\gamma_{\Gamma}^{N,2}$–to the single and double layer potentials corresponding to the wavenumber $k$ and a density $\varphi$ we define the four Helmholtz BIO:

[TABLE]

3.1 First kind (local) MTF

The main idea in MTF is to use Calderón projectors in each subdomain $\Omega_{j}$, $j=0,1,2$ and to use continuity of Dirichlet and Neumann traces. Indeed, by applying Green’s identities in a domain $\Omega_{j}$ we have that [10]

[TABLE]

For $j=0$, we first apply to the equation above the Dirichlet trace on $\partial\Omega_{0}$ to the domain $\Omega_{0}$ and obtain

[TABLE]

The continuity conditions of Dirichlet traces on $\partial\Omega_{0}$, that is

[TABLE]

can be expressed in the equivalent form

[TABLE]

where the operators $X_{01}$ and $X_{02}$ are defined as

[TABLE]

in terms of the restriction operators $\Pi_{ij}$ onto the interface $\Gamma_{ij}=\partial\Omega_{i}\cap\partial\Omega_{j}$–in the cases when $i\neq j$ and the domains $\Omega_{i}$ and $\Omega_{j}$ share an edge, and the extension-by-zero operators defined as

[TABLE]

for two functions $\psi_{01}$ and $\psi_{02}$ defined on $\Gamma_{01}$ and $\Gamma_{02}$ respectively. Using the reformulation (3.2) we can re-express equation (3.1) in the following form

[TABLE]

Remark 3.1

The expressions $u_{1}$ and $u_{2}$ in equation (3.3) should be understood as Dirichlet traces $u_{1}|_{\partial\Omega_{1}}$ and $u_{2}|_{\partial\Omega_{2}}$. We prefer the simpler notation as the meaning of those expressions should be grasped from the domain of definition of the operators that act upon them. Furthermore, all the restriction/extension operators can be defined for distributions [19].

Applying the normal derivative with respect to $\mathbf{n}_{0}$ to the same Green’s identities above we obtain

[TABLE]

The continuity conditions of Neumann traces on $\partial\Omega_{0}$

[TABLE]

can be expressed in the equivalent form

[TABLE]

and thus, equation (3.4) can be expressed as

[TABLE]

We define the Calderón block operators $\mathcal{C}_{j}$ satisfying

[TABLE]

and also introduce the following matrix operators:

[TABLE]

Denoting the Cauchy data on $\partial\Omega_{j}$ by

[TABLE]

and the incoming data $\mathbf{u}^{inc}$ on $\partial\Omega_{0}$ by

[TABLE]

we can re-write equations (3.3) and (3.6) in the form

[TABLE]

Similar techniques applied to $u_{1}$ in the domain $\Omega_{1}$ and respectively $u_{2}$ in the domain $\Omega_{2}$ lead to two additional equations in the vein of equation (3.7) leading to the local MTF [19]:

[TABLE]

where the additional matrix operators $\mathbb{X}_{j\ell}$ are similarly defined. The well-posedness of the local MTF in the multitrace space $\prod_{j=0}^{2}\left(H^{1/2}(\partial\Omega_{j})\times H^{-1/2}(\partial\Omega_{j})\right)$ was established in the same reference [19]. We note that owing to the presence of both weakly singular and hypersingular operators in the Calderón projectors $\mathcal{C}_{j},j=0,1,2$, MTF (3.8) are boundary integral equations of the first kind. As such, numerical solution of the MTF typically resort to preconditioning, an issue we return to in Section 6.

3.2 Second kind (global) MTF

An alternative boundary integral equation formulation of the transmission problem (2.1) can be derived using layer potentials defined on the skeleton $\Gamma:=\Gamma_{01}\cup\Gamma_{12}\cup\Gamma_{02}$ whose segments are oriented clockwise [8]. In this approach, the solutions $u_{j},j=0,1,2$ of the transmission problem (2.1) are sought in the form

[TABLE]

where $v$ and $p$ are densities defined on the skeleton $\Gamma$ and the double layer operators are defined with respect to exterior unit normals $\mathbf{n}$ with respect to the orientation chosen on $\Gamma$. Here we used the same convention that the index $j$ in equation (3.9) refers to the wavenumber $k_{j}$ for $j=0,1,2$. The enforcement of the continuity conditions of Dirichlet and scaled Neumann traces on the interfaces between subdomains leads to the following system of boundary integral equations, for all $\ell,j=0,1,2$ such that $\ell<j$:

[TABLE]

In what follows we refer to the formulation (3.2) by the acronym g-MTF. All of the boundary integral operators featured in equation (3.2) are defined on the skeleton $\Gamma$. The presence of (a) multiples of the densities $v$ and $p$ as well as (b) differences of hypersingular operators $N_{j}-N_{\ell}$ suggest referring to equations (3.2) as second kind boundary integral equations, hence the acronym CFIESK. We note that in contast to the MTF (3.8), the second kind equations (3.2) require one (scalar) unknown per interface $\Gamma_{j\ell}$. To the best of our knowledge, the well posedness of second kind boundary integral equations (3.2) is still an open question, see the discussion on this in [9]. Our numerical experiments indicate that these formulations appear to be well posed for a wide range of wavenumbers $k_{j},j=0,1,2$, thus corroborating the findings in [9].

4 Domain decomposition approach with classical Robin boundary conditions

DDM are natural candidates for numerical solution of transmission problems (2.1). A non-overlapping domain decomposition approach for the solution of equations (2.1) consists of solving subdomain problems with matching Robin boundary conditions on the common subdomain interfaces [11]. Indeed, this procedure amounts to computing the subdomain solutions:

[TABLE]

In equations (4.1), $\eta$ is assumed to be a positive number. The latter requirement is needed to ensure the well-posedness of the impedance boundary value Helmholtz problem in the exterior domain $\Omega_{0}$ [10].111In all the numerical examples in Section 6 we took $\eta=k_{0}$.

In order to describe the DDM method more concisely we introduce subdomain Robin-to-Robin (RtR) maps [16]. Given a subdomain $\Omega_{j}$ we define the RtR map $\mathcal{S}^{j}$ in the following manner:

[TABLE]

where $u_{j}$ is the solution of the following problem:

[TABLE]

In the case when $\Omega_{j}$ is the exterior domain $\Omega_{0}$, we further require in the definition of the RtR map $\mathcal{S}^{0}$ that $u_{0}$ is radiative at infinity. The DDM method computes the global Robin data $f=[f_{1}\ f_{2}\ f_{0}]^{\top}$ –the reader should observe the coordinate order– with

[TABLE]

as the solution of the following linear system that incorporates the subdomain RtR maps $\mathcal{S}^{j},j=0,1,2$, previously defined

[TABLE]

with right-hand side $g=[g_{1}\ g_{2}\ g_{0}]^{\top}$ wherein

[TABLE]

Remark 4.1

The domains $\Omega_{j}$, $j\neq 0$, can be further subdivided into smaller subdomains, in which case the DDM system (4.3) has to be augmented to incorporate the additional Robin data on the new interfaces. The size of the subdomains–in terms of wavelengths–should ideally be such that the computation/application of the corresponding RtR operators can be performed efficiently.

We note that the matrix $A$ in equation (4.3) is not stored in practice, and, due to its possibly large size, the DDM linear system (4.3) is typically solved in practice via iterative methods. Iterative solvers (e.g. Jacobi, GMRES) for the solution of DDM linear systems of the type described in equation (4.3) require large numbers of iterations, especially in the case of larger numbers of subdomains (e.g. see Table 5 in Section 6). This shortcoming can be attributed to the choice of classical Robin boundary conditions and the outflow/inflow of information from a subdomain to its neighboring subdomains associated with it. Ideally, the subdomain boundary conditions have to be chosen so that information flows out of the subdomain and no information is reflected back into the subdomain. This can be achieved if the term $i\eta$ in equations (4.3) is replaced by the adjacent subdomain Dirichlet-to-Neumann (DtN) operator restricted to the common interface—in this way the Jacobi scheme converges in precisely two iterations [27, 20]. Since DtN maps are not always well defined and expensive to calculate even when properly defined, easily computable approximations of DtN maps can be employed effectively to lead to faster convergence rates of GMRES solvers for DDM algorithms [4]. We present in Section 4.2 a non-overlapping DDM based on square root approximations of DtN operators.

In order to bypass the shortcomings of iterative solvers for DDM formulations (4.3), we resort to a direct solver that consists of application of recursive Schur complements. We begin by further splitting the subdomain data $f$ into interface components

[TABLE]

and by rewriting the linear sytem (4.3) in terms of the new unknown arrangement. To this end, we further decompose the RtR operators $\mathcal{S}^{1}$ and $\mathcal{S}^{2}$ in the block form

[TABLE]

where the block operators are defined informally as

[TABLE]

and

[TABLE]

We similalry split the operator $\mathcal{S}^{0}=[\mathcal{S}^{0}_{10}\quad\mathcal{S}^{0}_{20}]$. These block decompositions/splittings of the RtR maps $\mathcal{S}^{j},j=0,1,2$ allow us to write explicitly the DDM system (4.3) in the form:

[TABLE]

We note that the matrix in the new linear system (4.6) exhibits a block sparsity pattern that resembles that of matrices corresponding to two-dimensional finite difference operators that use five-point Laplaceans. We use a Schur complement approach for the solution of the linear system (4.6) that consists of elimination of the Robin data $(f_{12},f_{21})$ corresponding to the interior interface $\Gamma_{12}$ from the linear system (4.6). This procedure requires using the inverse of the matrix

[TABLE]

which can be computed explicitly

[TABLE]

Using the Schur complement in equation (4.6) we obtain the following linear system

[TABLE]

where $f_{0}^{int}=[f_{10}\quad f_{20}]^{\top}$. Interestingly, the operator $\mathcal{S}^{int}$ that features in equation (4.8) is itself an RtR operator defined as

[TABLE]

where $u$ is the solution of the following problem:

[TABLE]

and $u$ and $\varepsilon(x)^{-1}\partial_{n}u$ are continuous across $\Gamma_{12}$ where $\varepsilon(x)=\varepsilon_{1},\ x\in\Omega_{1},\ \varepsilon(x)=\varepsilon_{2},\ x\in\Omega_{2}$, and $n=n_{1}=-n_{2}$ on $\Gamma_{12}$. For obvious reasons, the procedure by which $\mathcal{S}^{int}$ was obtained above can be interpreted as the merging of subdomain RtR maps $\mathcal{S}^{j},j=1,2$.

Remark 4.2

In the case of five subdomain configuration depicted in Figure 2, the Schur complement elimination procedure proceeds by rearranging the subdomain Robin data in the order

[TABLE]

and then (1) first eliminating the unknown pairs $[f_{12},\ f_{21}]$ and $[f_{34},\ f_{43}]$ in the first stage, and (2) then subsequently eliminating the unknown pair $f_{top}:=[f_{23},\ f_{14}]$ and $f_{bottom}:=[f_{32},\ f_{41}]$. After performing the elimination stage we arrive at a sytem of the form (4.8) with $f_{0}^{int}=[f_{10}\ f_{20}\ f_{30}\ f_{40}]$ and the same interpretation of the operator $\mathcal{S}^{int}$ in the domain $\Omega_{1}\cup\Omega_{2}\cup\Omega_{3}\cup\Omega_{4}$.

The system (4.8) reduces in turn to the following equation involving the Robin unknown $f_{0}$ defined on the other interface $\partial\Omega_{0}$:

[TABLE]

Ideally, the linear system (4.10) should be solved by direct methods as well. Alternatively, if size of the exterior Robin unknown $f_{0}$ becomes too large for application of direct linear algebra solvers, Krylov subspace iterative solvers such as GMRES [29] can be employed for the solution of the reduced system (4.10). We note that the operator $I-\mathcal{S}^{int}\mathcal{S}^{0}$ is a pseudodifferential operator of order $-1$ (this can be rigorously established in the case when $k_{1}=k_{2}$), and thus equation (4.10) does not possess ideal spectral properties for fast convergence of iterative Krylov subspace solvers. Alternative interior/exterior coupling strategies that possess superior spectral properties is subject of ongoing investigations.

4.1 Analysis of the Schur complement DDM algorithm

As it can be seen from the description in Section 4, the DDM Schur complement procedure hinges on the invertibility of the operators $I-\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21}$, an issue which we examine in what follows. A more precise definition of the block operators in equations (4.4) and (4.5) can be given by considering the partial Robin-to-Robin maps:

[TABLE]

where

[TABLE]

Then we define

[TABLE]

and similarly for the other block operators.

The next results sheds light into the mapping properties of the operators $\mathcal{S}^{1}_{21,12}$ and $\mathcal{S}^{2}_{12,21}$:

Theorem 4.3

The following inequalities hold true

[TABLE]

Proof. This follow easily from the definition of the block operators $\mathcal{S}^{1}_{21,12}$ and $\mathcal{S}^{2}_{12,21}$ and the fact that $\mathcal{S}^{1}_{2}$ and $\mathcal{S}^{2}_{1}$ are isometries in $L^{2}(\partial\Omega_{1})$ and $L^{2}(\partial\Omega_{2})$ respectively. $\quad{}_{\blacksquare}$

We establish next the injectivity of the operator $I-\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21}$:

Theorem 4.4

The operator $I-\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21}:L^{2}(\Gamma_{12})\to L^{2}(\Gamma_{12})$ is injective.

Proof. Let $\psi\in\operatorname{Ker}(I-\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21})$, that is

[TABLE]

It is immediate to see that the operator $\mathcal{S}^{1}_{21,12}$ is invertible, and its inverse is defined as

[TABLE]

where

[TABLE]

with $u$ being the solution of the following boundary value problem:

[TABLE]

Define then $v_{2}$ to be the unique solution of

[TABLE]

and $v_{1}$ to be the unique solution of

[TABLE]

The relation $\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21}\psi=\psi$, which is equivalent to $\mathcal{S}^{2}_{12,21}\psi=(\mathcal{S}^{1}_{21,12})^{-1}\psi$ implies that

[TABLE]

The last equation together with the fact that

[TABLE]

imply that

[TABLE]

Defining then

[TABLE]

we see from equation (4.12) that $v$ is continuous across $\Gamma_{12}$, and, if we take into account that $n:=n_{1}=-n_{2}$ on $\Gamma_{12}$, so is $\varepsilon_{j}^{-1}\partial_{n}v$ across $\Gamma_{12}$. In addition, $v$ satisfies the following boundary value problem

[TABLE]

Taking into account the continuity properties of $v$ and its normal derivative across $\Gamma_{12}$, it follows immediately that

[TABLE]

which implies that $v=0$ on $\partial(\Omega_{1}\cup\Omega_{2})$. Given the zero impedance boundary condition on $\partial(\Omega_{1}\cup\Omega_{2})$, we also have that $\partial_{n}v=0$ on $\partial(\Omega_{1}\cup\Omega_{2})$, and hence $v=0$ in $\Omega_{1}\cup\Omega_{2}$. The last relation, in turn, implies that $\psi=0$ on $\Gamma_{12}$ which completes the proof. $\quad{}_{\blacksquare}$

The results in Theorem 4.3 and Theorem 4.4 do not guarantee the invertibility of the operator $I-\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21}$ in $L^{2}(\Gamma_{12})$. Indeed, the operators $\mathcal{S}^{1}_{21,12}$ and $\mathcal{S}^{2}_{12,21}$ are not strict contractions, nor are they compact in the aforementioned spaces. Nevertheless, one could attempt to replace the inverses of the operators $I-\mathcal{S}^{1}_{21,12}\mathcal{S}^{2}_{12,21}$–and its counterparts–needed in the DDM Schur complement algorithm by their truncated Neumann series expansions:

[TABLE]

Given that for small values of the parameter $N_{trunc}$ the evaluation of the expression in the right-hand side of equation (4.13) is less expensive than the full inversion, we explore in Section 6 the effectiveness of this approach.

4.2 DDM with generalized Robin boundary conditions

The rate of convergence of iterative Krylov subspace solvers of the DDM linear system (4.3) is largely determined by the choice of Robin boundary conditions therein. More effective Robin/impedance boundary conditions on the subdomain interfaces are known to improve the performance of iterative DDM solvers [30, 4, 34, 15]. These generalized Robin boundary conditions consist in replacing the classical $i\eta$ term by operators that approximate the Dirichlet-to-Neumann (DtN) operators of adjacent domains. For instance, the ideal Robin operator on the interface $\Gamma_{12}$ corresponding to the domain $\Omega_{1}$ consists of the operator $\varepsilon_{2}^{-1}Y^{2}|_{\Gamma_{21}}$, where $Y^{2}$ is the DtN operator corresponding to the domain $\Omega_{2}$ with zero Dirichlet boundary conditions on $\partial\Omega_{2}\setminus\Gamma_{21}$. With this very choice, the ensuing DDM algorithm converges in precisely two iterations [27], at least in the case when $\Omega_{j}$ are half planes. Similarly, the ideal Robin operator on the interface $\Gamma_{10}$ corresponding to the domain $\Omega_{1}$ can be shown to consist of the operator $\varepsilon_{0}^{-1}Y^{0}|_{\Gamma_{01}}$. Although DtN operators are not always defined for interior subdomains, they are always well defined in the exterior domain $\Omega_{0}$. However, even when they are properly defined DtN maps are non-local operators whose computation can be expensive. Their computation, whenever possible, can be obtained via boundary integral operators. For instance, using Green’s identities in the domain $\Omega_{2}$ and taking into consideration the null Dirichlet boundary conditions on $\partial\Omega_{2}\setminus\Gamma_{21}$, we obtain the expression:

[TABLE]

which upon application of Dirichlet traces, leads to the identity

[TABLE]

The invertibility of the operator $S_{\partial\Omega_{2},2}$ in the equation above, and hence the well-posedness of the DtN operator $Y^{2}$, can be guaranteed provided the subdomain $\Omega_{2}$ is small enough, typically less than one wavelength across. A simple solution that would allow one to consider subdomains of any size is to consider DtN operators $Y^{2,c}$ corresponding to complexified wavenumbers $k_{2}+i\sigma_{2},\ \sigma_{2}>0$ instead of the operators $Y^{2}$. Using these operators, we can define a transmission operator on the interface $\partial\Omega_{2}$ in the form

[TABLE]

and similar transmission operators on the interfaces $\partial\Omega_{2}$ and $\partial\Omega_{0}$ respectively. We then match DtN Robin boundary conditions (DtNR) on the subdomain interfaces in the form

[TABLE]

Remark 4.5

We point out several issues related to the DDM formulation that uses transmission operators defined in equation (4.15). The primary issue related to generalized Robin boundary conditions matching is whether these are equivalent to the continuity of the Dirichlet and (scaled) Neumann traces on subdomain interfaces. It is a straightforward matter to see that the aforementioned equivalence is guaranteed as long as the operators $\varepsilon_{2}^{-1}\Pi_{21}Y^{2,c}+\varepsilon_{1}^{-1}\Pi_{12}Y^{1,c}$ and their conterparts are invertible in appropriate functional spaces of functions/distributions on $\Gamma_{12}$. These functional spaces are related to the mapping property $\Pi_{21}Y^{2,c},\Pi_{12}Y^{1,c}:\widetilde{H}^{1/2}(\Gamma_{12})\to H^{-1/2}(\Gamma_{12})$ (we used the same notation and definitions for the Sobolev spaces on the interfaces $\Gamma_{12}$ as in reference [26]). Following similar arguments to those in the proof of Theorem 4.4, we can establish the invertibility of the operator $\varepsilon_{2}^{-1}\Pi_{21}Y^{2,c}+\varepsilon_{1}^{-1}\Pi_{12}Y^{1,c}$. The second issue is related to the mapping property of the operators $\mathcal{T}_{1}^{DtN}$ defined in equation (4.15) and the well-posedness of the Helmholtz equation in the domain $\Omega_{1}$ with boundary conditions that feature the generalized Robin operator of the left-hand side of equation (4.16). Clearly, the straightforward extension by zero of the operator $\Pi_{21}Y^{2,c}:\widetilde{H}^{1/2}(\Gamma_{12})\to H^{-1/2}(\Gamma_{12})$ (e.g. the operator $X_{12}Y^{2,c}$) does not give rise to an operator that maps $H^{1/2}(\partial\Omega_{1})$ to $H^{-1/2}(\partial\Omega_{1})$, and thus, the operators $\mathcal{T}_{1}^{DtN}$ defined in equation (4.15) do not enjoy the latter mapping property. We discuss in what follows a localization/extension procedure via smooth cut-off functions that delivers transmission operator with desired mapping properties.

Besides the issues raised in Remark 4.5, the DtN operators for the domains $\Omega_{j},0\leq j$, even when properly defined, are expensive to compute and to apply in equations (4.16). We therefore rely on suitable approximations of DtN operators that can be readily computed and at the same time lead to well posed generalized Robin boundary conditions in the domains $\Omega_{j},\ 0\leq j$. These approximations start from an exact representation of the DtN operators in each domain. For instance, an application of Neumann traces to the Green’s identities in the domain $\Omega_{2}$ leads to

[TABLE]

and thus, whenever admissible,

[TABLE]

An explicit approximation of the DtN operator $Y^{2}$ can be achieved by keeping the leading order term in the equation above in the sense of pseudodifferential operators

[TABLE]

The restriction of the operators $N_{\partial\Omega_{2},2}$ to the interface $\Gamma_{12}$–which is an open arc–has to be performed with care [31, 33]. In addition, the generalized Robin operators should lead to well posed problems in the domain $\Omega_{1}$. One possibility to blend the approximations of $Y^{2}|_{\Gamma_{21}}$ and $Y^{0}|_{\Gamma_{20}}$ that use hypersingular operators is given by [33]:

[TABLE]

where $\chi_{01}$ and $\chi_{12}$ are smooth cutoff functions whose supports lie in $\Gamma_{10}$ and $\Gamma_{12}$, respectively. Given that hypersingular operators are expensive to compute, we can replace the hypersingular operators in equation (4.17) by principal symbol Fourier multiplier operators. The latter principal symbols are defined as

[TABLE]

where the square root branches are chosen such that the imaginary parts of the principal symbols are positive. The principal symbol Fourier multipliers are defined in the Fourier space $TM(\partial\Omega_{1})$ [1] as

[TABLE]

for a density $\varphi_{1}$ defined on $\partial\Omega_{1}$. We define accordingly

[TABLE]

and we use the operators in equation (4.19) as generalized Robin/impedance operators. Thus, the classical Robin boundary conditions on $\partial\Omega_{1}$ are replaced by the following Generalized Square Root Robin boundary conditions (GSqR)

[TABLE]

We note that the operators $PS(\mathcal{T}_{1}):H^{1/2}(\partial\Omega_{1})\to H^{-1/2}(\partial\Omega_{1})$ are coercive, and thus the Helmholtz equations in the domain $\Omega_{1}$ with generalized Robin boundary conditions defined by the operators in the left-hand side of equation (4.20) are well posed [33]. Similar generalized impedance operators can be defined for the domains $\Omega_{0}$ and $\Omega_{2}$ and then incorporated in a DDM algorithm that computes the generalized Robin data

[TABLE]

by making use of suitably defined generalized RtR maps $\mathcal{S}^{g,j}$.

Remark 4.6

While the new blended operators $PS(\mathcal{T}^{j}),j=0,1,2$ give rise to well-posed Helmholtz equations in each subdomain, the issue of the invertibility of the operator $(PS(\mathcal{T}^{1})+PS(\mathcal{T}^{2}))|_{\Gamma_{12}}$–and hence the equivalence between conditions (4.20) and the continuity of the Dirichlet and Neumann traces on subdomain interfaces–is an open question. In addition, the generalized RtR maps are no longer unitary and thus the proof of solvability of the DDM linear system with GSqR (4.20) also remains open.

5 High-order Nyström discretizations of the MTF equation (3.8), g-MTF equation (3.2) and the DDM formulation (4.3)

We use Nyström discretizations of the MTF equation (3.8), g-MTF equation (3.2) and the DDM formulation (4.3) that relies on discretizations of the four BIO in the Calderón calculus. The latter, in turn, encompass: (a) use of graded meshes based on sigmoid transforms that cluster polynomially discretization points toward corners; (b) splitting of the kernels of the parametrized versions of the boundary integral operators that feature in the MTF equation (3.8) and the g-MTF equation (3.2) into sums of regular quantities and products of periodized logarithms and regular quantities; (c) trigonometric interpolation of the densities of the boundary integral operators; and, (d) analytical expressions for the integrals of products of periodic singular and weakly singular kernels and Fourier harmonics. These discretizations were introduced in [12] where this methodology is presented in full detail. The main idea of our Nyström discretization is to incorporate sigmoid transforms [23] in the parametrization of a closed Lipschitz curve $\Gamma$ and then split the kernels of the Helmholtz BIO into smooth and singular components. Using graded meshes that avoid corner points and classical singular quadratures of Kusmaul and Martensen [24, 25], we employ the Nyström discretization presented in [12] to produce high-order approximations of the BIO that enter the Calderón projectors. We note that the role of weighted traces is to increase the regularity of the densities to be integrated, regularity that can be tuned by simply increasing the polynomial order in the sigmoid transforms. In what follows, we provide details relevant to each of the three formulations considered.

For each of the subdomains $\Omega_{j}$, $j=0,1,2$, we employ graded meshes $\mathbf{x}^{j}_{m},m=0,\ldots,N_{j}-1$ on $\partial\Omega_{j}$ with the same polynomial degree of the sigmoid transform. All these meshes are subsequently shifted by the same amount–basically half a discretization step-size–so that none of the grid points corresponds to a triple/multiple junction or a corner point. We note that on a common interface between two subdomains $\Omega_{j}$ and $\Omega_{\ell}$ the grid points corresponding to the mesh in each subdomain coincide.

Nyström discretization of the MTF equation (3.8). Following the prescriptions in [12] for the discretization of transmission boundary integral equations for one subdomain, i.e. no triple/multiple junctions, we define weighted Neumann traces on the boundary of each subdomain $\Omega_{j}$

[TABLE]

where $\mathbf{x}_{j}$ is a $[0,2\pi]$ parametrization of the boundary $\partial\Omega_{j}$ that incorporates sigmoid transforms on each of the smooth arcs of $\partial\Omega_{j}$. This procedure requires introducing weighted parametrized versions of the four BIO of the Calderón calculus; these were discussed in detail in the same reference [12].

Nyström discretization of the g-MTF equation (3.2). The additional challenge in the Nyström discretization of the g-MTF consists of the fact that the skeleton $\Gamma$ is no longer a closed curve–recall that the Nyström discretizations that rely on the singular quadratures of Kusmaul and Martensen [24, 25] require that the domain of integration is a closed curve and thus all the integrands to be dealt with are periodic. Nevertheless, introducing weighted quantities

[TABLE]

we obtain a weighted g-MTF whose unknowns are $(v^{w},p^{w})$ by multiplying both sides of equations (3.2) by the weight $|\mathbf{x}^{\prime}|$. We note that the weight $|\mathbf{x}^{\prime}|$ is related to a parametrization of the skeleton $\Gamma$ that is easily produced from individual parameterizations of the domains $\partial\Omega_{j},j=0,1,2$. We note that the weighted unknowns $(v^{w},p^{w})$ vanish algebraically at the triple/multiple junction points of the skeleton $\Gamma$. This simple fact allows us to evaluate easily weighted boundary integral operators defined on $\Gamma$. Indeed, the evaluation of a parametrized weighted integral operator of the type

[TABLE]

can be simply performed as

[TABLE]

where

[TABLE]

Each of the two integrals in equation (5.1) is performed on a closed curve and involves $2\pi$ periodic and regular functions, and thus can be computed with high-order accuracy by the Nyström method depicted above. All of the weighted integral operators that enter the weighted g-MTF can be effected in a similar manner.

Discretization of the DDM equation (4.3) based on Nyström discretizations of the RtR maps. The RtR maps $\mathcal{S}^{j}$ can be computed via well-conditioned BIO and subsequently discretized following the Nyström discretizations introduced in [33]. In a nutshell, the method introduced in [33] is a direct regularized boundary integral equation formulation of interior/exterior Helmholtz equations with Robin/impedance boundary conditions of the type employed by the DDM formulation (4.3). Specifically, if $u_{j}$ is the solution of the following Helmholtz problem with Robin boundary conditions, augmented by $u_{0}$ is radiative in the exterior domain $\Omega_{0}$,

[TABLE]

then $u_{j}|_{\partial\Omega_{j}}$ can be shown to satisfy the following well-conditioned regularized boundary integral equation [33]:

[TABLE]

where $Z_{j}=i\eta\varepsilon_{j}$. Thus the RtR operator $\mathcal{S}^{j}$ can be expressed as

[TABLE]

The generalized Robin conditions (4.15) and (4.20) are amenable to the same treatment provided in [33], the only difference being that the scalar impedance $Z_{j}$ in equations (5) and (5.3) is replaced by impedance operators defined in equations (4.16) and (4.19) respectively. In order to avoid complications related to singularities at junction/cross points, we replace in the DDM algorithm the RtR maps by weighted parametrized counterparts

[TABLE]

Collocated discretizations of the latter weighted RtR maps can be easily computed through a simple modification of the methodology introduced in [33] and recounted above. However, the representation of RtR maps in terms of BIO requires use of inverses of the operators $\mathcal{A}_{j}$ cf. equation (5.3). In order for the DDM algorithm to be efficient, the electromagnetic/acoustic size of subdomains $\Omega_{j}$ should be amenable to application of direct linear algebra solvers for calculations of the inverses of the collocation of the matrices $\mathcal{A}_{j}$. The discretization of the weighted RtR maps corresponding to each domain $\partial\Omega_{j}$ is thus constructed as $N_{j}\times N_{j}$ collocation matrices $\mathcal{S}^{j}_{N}$. Alternatively, the application of the collocated versions of the operators $\mathcal{S}^{j}$ can be effected via Krylov subspace iterative solutions of equations (5) for each domain $\Omega_{j}$. However, the numbers of iterations needed for solution of equations (5) does grow with the electromagnetic/acoustic size of the interior bounded domains $\Omega_{j}$, $j\neq 0$ [33]. Therefore, it is advisable to further subdivide the interior domains $\Omega_{j}$ when their size is deemed too large for efficient computations of their corresponding RtR maps. In the case of the exterior/unbounded domain $\Omega_{0}$, the numbers of iterations required for the solution of equations (5) are small regardless of the frequency $\omega$ [33], and thus the application of $\mathcal{S}^{0}_{N}$ can be performed effciently by iterative methods.

Once the collocation matrices $\mathcal{S}^{j}_{N}$ are constructed, the discretization of the interface subdomain RtR maps (e.g. the maps $\mathcal{S}^{j}_{j\ell,\ell j}$ and all the other ones defined in equation (4.4)) is straightforward since it simply amounts to extracting suitable blocks from the matrices $\mathcal{S}^{j}_{N}$. Indeed, the discretization of the operators $\mathcal{S}^{j}_{j\ell,\ell j}$ consists of extracting from the collocation matrix $\mathcal{S}^{j}_{N}$ the block that corresponds to self-interactions of the grid points on the interface/edge $\partial\Omega_{j}\cap\partial\Omega_{\ell}$. This is possible since none of these mesh points $\mathbf{x}^{j}_{m}$ corresponds to a corner of $B_{j}$. We also mention that all the matrix inversions needed in the Schur complement algorithm (cf. formula (4.7)) are performed by either direct linear algebra methods: the size of these matrices does not exceed the number of unknowns need to discretize $\partial\Omega_{0}$ or by Neumann series (cf. equation (4.13)). In case when the large domains $\Omega_{j}$, $j\neq 0$, are further subdived into smaller non-overlapping subdomains $\Omega_{j}=\cup_{n=1}^{J}\Omega_{j_{n}}$ for which RtR maps $\mathcal{S}^{j_{n}}_{N_{j_{n}}}$ can be computed via direct linear algebra methods. The RtR maps $\mathcal{S}^{j}_{N}$ themselves can be computed by merging the maps $\mathcal{S}^{j_{n}}_{N_{j_{n}}}$. However, as pointed in Table 5 in Section 6, the numbers of iterations needed by an iterative solver of the classical DDM system (4.3) as well as the DtNR and GSqR DDM system do grow with the number of additional subdomains.

The comments above are also applicable to the DtNR and GSqR DDM algorithms: the size of the interior subdomains has to be amenable to direct solvers computations of the generalized RtR maps. The DtNR DDM algorithm require that the subdomain DtN maps are precomputed, which is also done via Nyström solvers of the boundary integral formulation (4.14). The discretization of the operators $\mathcal{T}_{j}^{DtN}$defined in equations (4.15) is obtained by simply extracting suitable blocks from the collocation matrices corresponding to discretized DtN maps. These blocks are subsequently used in the computation of collocated generalized RtR maps. We note that given that the Fourier multiplier operators defined in equations (4.19) can be efficiently discretized via FFTs [12], the cost of computing subdomain generalized RtR maps is virtually identical to that of computing subdomain classical RtR map. This is in contrast to volumetric solvers for local domain problems (e.g. finite difference/finite element) where Padé approximations need to be employed [34] for approximating the square root operators in equations (4.19). Both DtNR and GSqR DDM algorithms require significantly fewer Krylov subspace iterations than the classical DDM algorithm, and the reduction in number of iterations is more striking in the high-contrast cases.

6 Numerical results

We present in what follows numerical experiments with various formulations of the transmission problems (2.1) considered in this text. We use as test configurations with piece-wise constant dielectric permittivities $\varepsilon_{j}$ in geometries depicted in Figures 1, in which the circular domain has radius one. All of the formulations considered were discretized following the prescriprion in Section 5. We used discretizations that place the same number of discretization points on each interface–matching boundary meshes, and sigmoid transforms of polynomial degree 3 for the MTF equation (3.8) and the DDM formulation (4.3) with $\eta=k_{0}$, and respectively degree 4 for the g-MTF (3.2). We chose $\sigma_{j}=0.1$ in the definition of the operators (4.15) incorporated in the DtNR and, respectively, $\sigma_{j}=k_{j}^{1/3}$ in the definition of the operators (4.19) that are incorporated in GSqR. We mention that the more robust choice of operators in equation (4.15) gives rise to DDM whose iterative performance is identical to that of DDM that use exact DtN maps. In all the numerical experiments we present numbers of GMRES iterations for various solvers to reach a residual of $10^{-4}$ and present errors in the far-field for $1,024$ equi-spaced far-field directions. The far-field errors were computed using fine discretizations of the MTF solver. As previously mentioned, the discretizations of both g-MTF and the DDM formulation use the same number of unknowns, whereas that of MTF uses twice as many unknowns.

On account of MTF formulation (3.8) being an integral formulation of the first kind, the number of iterations needed by Krylov subspace iterative solvers (e.g. GMRES) of its discretization can be considerable, especially for high-contrast configurations at high-frequencies. Fortunately, simple and quite effective preconditioners of the MTF formulations are readily available. Indeed, a simple preconditioner of the MTF fomulation (3.8) consists of the Calderón diagonal preconditioner [19]:

[TABLE]

Another possibility is to eliminate the unknown $\mathbf{u}_{0}$ from the MTF formulation (3.8) using Schur complements. Indeed, using the fact that $\mathcal{C}_{0}$ is a projector, we get that $\mathbf{u}_{0}$ can be explictly eliminated from the MTF system (3.8) via the equation

[TABLE]

leading to a reduced MTF formulation:

[TABLE]

Formulation (6.2), in turn, can be further preconditioned by the operator $diag(\mathcal{C}_{1},\mathcal{C}_{2})$, a formulation for which we report the numbers of GMRES iterations in this section.

We start the presentation of numerical results by illustrations in Tables 1 and 2 of the performance of the three types of solvers dicussed in this text in the case of multiple junction configurations with very high-contrast materials at low incident frequencies $\omega$. We continued in Tables 3 and 4 with results concerning multiple junction configurations at high frequencies. Finally, we conclude in Table 5 with an illustration of the behavior of the Nyström solvers of the formulations considered as the size of the discretization is increased. The accuracy achieved by the GSqR DDM solver is at the same level to that of the DDM solver in all of the cases considered, and for this reason and in order to better streamline the information presented we chose not to report it.

Discussion of the performance of various formulations. We begin by noting that the computational cost of setting up the MTF equation (3.8) and the g-MTF (3.2) is comparable, as the building blocks of these formulations are the same BIO defined on the boundaries of the subdomains. The behavior of iterative solvers for the solution of Nyström discretizations of the two formulations is quite different, given that MTF is a first kind formulation that requires twice as many unknowns than the second kind formulation CFIESK. Calderón preconditioning strategies that lead to the MTF formulation (6.1) and Schur complement strategies combined with Calderón preconditioning strategies that lead to the MTF formulation (6.2) are effective in reducing the numbers of GMRES iterations and do not require significant additional computational costs. In contrast, the g-MTF (3.2), whose motivation was canceling the adverse effects of hypersingular operators, is not easily amenable to preconditioning. However, even with preconditioning, both formulations require large numbers of GMRES iterations, especially in the case of high-contrast multi-domain configurations at high-frequencies.

As it can be seen from the results presented in this section, solvers based on DDM formulations (4.3) and the DDM formulation that incorporate either the DtNR or GSqR give rise to smaller numbers of GMRES iterations than those based on the MTF and g-MTF. We note that effecting the matrix-vector product associated with the DDM equation (4.3), requires in turn application of each of the RtR maps $\mathcal{S}^{j}$. The latter can be done iteratively via solutions of Robin/impedance boundary values problems in each corresponding domain $\Omega_{j}$ via equations (5). We present in Table 6 the numbers of GMRES iterations required for the solution of equations (5) in each subdomain $\Omega_{j}$ when using a residual of $10^{-4}$. We observe that these numbers grow with the frequency $\omega$, especially in the case of the interior domains $\Omega_{1}$ and $\Omega_{2}$. Further subdivisions of the interior domains bring down these numbers, but only modestly so, while increasing considerably the numbers of GMRES iterations needed for the solution of the global DDM system (4.3).

These findings corroborate those reported in the DDM literature [30] and suggest that DDM algorithm are efficient when the size of subdomains is amenable to use of direct solvers for computation of the action of subdomains RtR maps. The Nyström discretizations presented in Section 5 constitute very efficient means to discretize both classical and generalized RtR operators. Our DDM algorithms thus proceed in two stages: (1) a precomputation/offline stage where all the subdomain RtR maps are precomputed, and (2) an online stage where the DDM linear system is solved iteratively in a matrix-free fashion. We note that the use of GSqR (4.19) in the DDM algorithm gives rise to important savings in numbers of iterations with virtually no additional cost. The numbers of iterations is further decreased in the DtNR DDM algorithm, yet at the expense of adding one more offline stage (0) where the subdomain DtN maps are precomputed. However, the number of subdomains adversely affects the rate of convergence of iterative solvers of the DDM system (4.3), regardless of the type of Robin boundary conditions used. These findings constitute somewhat of a surprise especially in the case when DtNR boundary conditions (4.16) are used: one would expect that the incorporation of (almost) exact DtN should result in numbers of DDM iterations that are insensitive to the frequency and the number of subdomains. A similar situation was analyzed in the one dimensional case with constant wavenumber in [34]). The use of double sweep preconditioners in slab-type DDM was shown in [34] to scale with the frequency and the number of subdomains in the absence of sharp interface sof material discontinuity. The incorporation of double sweep preconditioners for the GSqR DDM is currently underway.

The Schur complement DDM algorithm does require that the collocation matrices $\mathcal{S}^{j}_{N}$ be precomputed for all the interior subdomains, i.e. $j\neq 0$. This can be done by either inverting matrices of sizes $N_{j}\times N_{j}$, $j\neq 0$ at a cost $N_{j}^{3}$, or, alternatively, by further subdividing the domain $\Omega_{j}$ into smaller subdomains –for which the corresponding collocated RtR maps can be computed effciently with direct methods– and then applying the merging procedure to the RtR maps of the smaller subdomains. The latter algorithm can also be computationally more advatangeous than the direct approach, especially in the case of identical subdomains. The Schur complement DDM algorithm requires at each stage inversion of matrices (cf. formula (4.7)) whose size depends on the number of unknowns on the common interfaces between subdomains. The matrix inversions required by the merging procedure can be done using direct linear algebra method or using Neumann series (4.13). If the interior domain $\cup_{1\leq j}\Omega_{j}$ is decomposed into a collection $P\times P$ subdomains, and $n$ discretization points are used per each subdomain interface, then the cost of hierarchically merging of the subdomain RtR maps to deliver the interior RtR map $\mathcal{S}^{int}$ is $\mathcal{O}(P^{3}n^{3})$. In contrast, direct solvers for MTF would require computational costs of the order $\mathcal{O}(2^{9}P^{6}n^{3})$, whereas direct solvers for g-MTF would require costs of the order $\mathcal{O}(2^{6}P^{6}n^{3})$. If the aforementioned matrix inversions are bypassed via truncated Neumann series, the number of terms needed in the Neumann series is quite small in the case of low contrast/low frequencies, but does grow in the case of very high-contrast/ high frequencies.

Finally, after the merging of the RtR maps of interior domains $\mathcal{S}^{int}$ has been computed via hierarchical Schur complements, the equation (4.10) is solved either via direct linear algebra solvers or iteratively. If the latter approach is used, the exterior RtR map $\mathcal{S}^{0}$ needs not be precomputed. Furthermore, the computation of a matrix vector product associated with the exterior RtR map $\mathcal{S}_{0}$ can be performed efficiently via iterative solutions of equation (5) according to the results in Table 6. As it can be seen from the results presented in this section, the number of GMRES iterations needed to solve equation (4.10) is quite reasonable throughout the range of problems considered, and does not appear to grow dramatically with the frequency. Nevertheless, equation (4.10) shares the same features with integral equations of the first kind, alternative robust interior/exterior coupling strategies are desirable. One possibility currently under investigation is to modify the generalized Robin operators (4.19) so that the same complex wavenumber is used on all the interior interfaces. We conclude by mentioning that, just like the classical iterative DDM algorithm, the Schur complement DDM algorithm is highly parallelizable.

7 Conclusions

We presented a comparison between the performance of Krylov subspace iterative solvers based on Multiple Trace Formulations (MTF) and Domain Decomposition Methods (DDM) for two dimensional frequency domain scattering problems involving bounded composite materials,i.e. piece-wise constant material parameters. We investigated DDM based both on classical and generalized/DtN Robin boundary conditions. The former DDM version gives rise to linear systems which, although not particularly well-suited for iterative solvers, can be solved efficiently via Schur complements. On the other hand, the generalized DDM version that incorporates square root approximations of DtN maps gives rise to relatively small numbers of iterations. Extensions to three-dimensional configurations are currently underway.

Acknowledgments

Catalin Turc gratefully acknowledges support from NSF through contracts DMS-1312169 and DMS-1614270. Carlos Jerez-Hanckes thanks partial support from Conicyt Anillo ACT1417.

8 Appendix

In this section we consider DDM for the Helmholtz equation in one dimension. More precisley, we consider the Helmholtz equation

[TABLE]

where the wavenumber $k(x)$ is a piecewise constant function, that is

[TABLE]

and $u$ and $u^{\prime}$ are continuous at $a_{j},j=0,\ldots,N+1$. We note that we do not require that the wavenumbers $k_{j}$ be necessarily different on adjacent intervals. The classical DDM formulation of the Helmholtz equation above can be written in the form

[TABLE]

for all $1\leq j\leq N$ together with the end-interval equations

[TABLE]

and

[TABLE]

To each of these Robin problems we associate RtR maps. First we define for $1\leq j\leq N$ the following matrices

[TABLE]

then the following complex scalars

[TABLE]

and

[TABLE]

Denoting $h_{j}=a_{j+1}-a_{j}$, it is a straightforward matter to compute

[TABLE]

and

[TABLE]

for $1\leq j\leq N$. We also get

[TABLE]

and

[TABLE]

Ordering the data $f=[f_{01}\ f_{10}\ f_{12}\ \ldots f_{N+1,N}]^{\top}$, then the classical DDM can be written in the form $(I+A)f=g$ where the matrix $I+A$ is given in explicit form

[TABLE]

We present in Figure 4 the spectral properties of the matrix $I+A$ defined in equation (8.2) for a case of piecewise constant wavenumber that takes four values in the interval $(0,1)$, and a total of 300 subintervals. The spectral properties of the ensuing DDM are associated with poor behavior of GMRES iterative solvers: the eigenvalues are distributed almost uniformly on a circle of radius close to one centered at $(1,0)$.

In the case of DtN DDM algorithm we make use of the following DtN maps, assumed to be properly defined:

[TABLE]

where $v_{j}$ is the solution of the following problem

[TABLE]

and

[TABLE]

where $w_{j}$ is the solution of the following problem

[TABLE]

It can be easily shown that

[TABLE]

The DtN DDM formulation of the Helmholtz equation above can be written in the form

[TABLE]

for all $1\leq j\leq N$ together with corresponding end-interval equations. Corresponding RtR DtN maps/matrices can be defined and their entries are given by

[TABLE]

and

[TABLE]

where

[TABLE]

With these notation the DtN DDM can be written in the form $(I+A^{dtn})f^{dtn}=g^{dtn}$ where the matrix $I+A^{dtn}$ is given in explicit form

[TABLE]

In what follows we present a strategy to construct an effective preconditioner for the linear systems associated with the classical DDM formulations. The main idea is to replace the entries in the matrix $I+A^{dtn}$ corresponding to $\mathcal{S}^{dtn,j}_{11}$ and $\mathcal{S}^{dtn,j}_{22}$ by zero. At this stage we find more intuitive to refer to $\mathcal{S}^{dtn,j}_{12}$ as to $\mathcal{S}^{dtn,j}_{b}$ (the subscript stands for backward, consistent with the direction in which the information propagates) and to $\mathcal{S}^{dtn,j}_{21}$ as to $\mathcal{S}^{dtn,j}_{f}$ (the subscript stands for forward). We consider thus the following matrix

[TABLE]

The matrices $I+\widetilde{A}^{dtn}$ corresponding to the same experiment described in Figure 4 have only one eigenvalue $\lambda=1$ with algebraic multiplicity $2(N+1)$ (this is the number of unknown in the DDM) and geometric multiplicity 2, that is it has only two linearly independent eigenvectors, which turn out to be the first and the last canonical vectors in $\mathbb{R}^{2(N+1)}$. This situation was already pointed out in [34] in the case of constant wavenumber. The key property that facillitates the use of the matrix $I+\widetilde{A}^{dtn}$ as the basis of a DDM preconditioner is the fact that the inverse of this matrix can be computed explicitly, and its expression does not involve algebraic inverses.