A Semi Matrix-Free Twogrid Preconditioner for the Helmholtz Equation with Near Optimal Shifts

Abstract

Due to its significance in terms of wave phenomena a considerable effort has been put into the design of preconditioners for the Helmholtz equation. One option to design a preconditioner is to apply a multigrid method on a shifted operator. In such an approach, the wavenumber is shifted by some imaginary value. This step is motivated by the observation that the shifted problem can be more efficiently handled by iterative solvers when compared to the standard Helmholtz equation. However, up to now, it is not obvious what the best strategy for the choice of the shift parameter is. It is well known that a good shift parameter depends sensitively on the wavenumber and the discretization parameters such as the order and the mesh size. Therefore, we study the choice of a near optimal complex shift such that a flexible generalized minimal residual (FGMRES) solver converges with fewer iterations. Our goal is to provide a map which returns the near optimal shift for the preconditioner depending on the wavenumber and the mesh size. In order to compute this map, a data driven approach is considered: We first generate many samples, and in a second step, we perform a nonlinear regression on this data. With this representative map, the near optimal shift can be obtained by a simple evaluation. Our preconditioner is based on a twogrid V-cycle applied to the shifted problem, allowing us to implement a semi matrix-free method in which only the coarse grid and boundary matrices need to be stored in memory. On the fine grid, only the action of the matrix applied to a vector is computed, without assembling the global matrix. This enables efficient use of computing resources and allows problems to be solved at scales that were previously limited by the available memory. The performance of our preconditioned FGMRES solver is illustrated by several benchmark problems with heterogeneous wavenumbers in two and three space dimensions.

A Local Fourier Analysis for the Twogrid Operator

The convergence analysis by means of LFA is based on the assumption that the error function $e_h^l$ on the fine grid $\textbf{G}_h$ can be represented by a linear combination of Fourier modes [41]:

$$\begin{aligned} \varphi _h( {\theta },\textbf{x} ) = e^{\textrm{i}{\theta } \cdot \textbf{x}/h} = e^{\textrm{i}\theta _1 x_1/h} e^{\textrm{i}\theta _2 x_2/h},\;\textbf{x} \in \textbf{G}_h. \end{aligned}$$

$\theta = ( \theta _1,\theta _2 ) \in \mathbb {R}^2$ denote the Fourier frequencies, which may be restricted to the domain $\varTheta = \left( -\pi ,\pi \right] ^2 \subset \mathbb {R}^2$, see [12], as a consequence of the fact that

$$\begin{aligned} \varphi _h( {\theta } + 2\pi ,\textbf{x} ) = \varphi _h( {\theta },\textbf{x} ). \end{aligned}$$

Thus the space of Fourier modes is given by:

$$\begin{aligned} E_h = \text {span}\left\{ \left. \varphi _h( {\theta },\textbf{x} ) = e^{\textrm{i}{\theta } \cdot \textbf{x}/h} \right| \textbf{x} \in \textbf{G}_h,\; {\theta } \in \varTheta \right\} . \end{aligned}$$

According to [41], the Fourier modes are the formal eigenfunctions of an operator

$$\begin{aligned} H \in \left\{ T_h^{2h},\;I_h,\;I_h^{2h},\;I_{2h}^h,\;L_{h},\;L_{2h},\;S_h \right\} . \end{aligned}$$

This means that:

$$\begin{aligned} H \varphi _h( {\theta },\textbf{x} ) = \tilde{H}( {\theta } ) \varphi _h( {\theta },\textbf{x} ), \end{aligned}$$

where $\tilde{H}( {\theta } )$ is called the Fourier symbol of the operator H. In a next step, we decompose the frequency space into a low frequency space $T^{low} = \left( -\frac{\pi }{2},\frac{\pi }{2}\right] ^2$ and a high frequency space $T^{high} = \varTheta \setminus T^{low}$. For each low frequency ${\theta } \in T^{low} = \left( -\frac{\pi }{2}, \frac{\pi }{2} \right] ^2,$ we define the following four frequencies:

$$\begin{aligned} {\theta }^{(0,0)} = ( \theta _1,\theta _2 ),\; {\theta }^{(1,1)} = ( \overline{\theta }_1,\overline{\theta }_2 ),\; {\theta }^{(1,0)} = ( \overline{\theta }_1,\theta _2 ),\; {\theta }^{(0,1)} = ( \theta _1,\overline{\theta }_2 ), \end{aligned}$$

where

$$\begin{aligned} \overline{{\theta }}_i = {\left\{ \begin{array}{ll} {\theta }_i + \pi &{}\quad \text { if } \theta _i<0,\\ {\theta }_i - \pi &{}\quad \text { if } \theta _i \ge 0. \end{array}\right. } \end{aligned}$$

Obviously, the four frequencies fulfil

$$\begin{aligned} \varphi _h( {\theta }^\alpha ,\textbf{x} ) = \varphi _{2h}( 2{\theta }^{(0,0)},\textbf{x} ),\; \textbf{x} \in \textbf{G}_{2h}, \;\mathbf {\alpha } \in I = \left\{ (0,0),\;(1,1),\;(1,0),\;(0,1)\right\} . \end{aligned}$$

Fourier modes with this relationship are called harmonic to each other. Considering a given low frequency ${\theta } = {\theta }^{(0,0)}$, we define a four dimensional space of harmonics by

$$\begin{aligned} E_h^{{\theta }} = \text {span}\left\{ \left. \varphi ( {\theta }^{\alpha }, \cdot ) \right| \mathbf {\alpha } \in I \right\} . \end{aligned}$$

An important feature of these spaces is that they are invariant under the twogrid operator $T_h^{2h}$. Applying $T_h^{2h}$ to an arbitrary element $\varPsi \in E_h^{{\theta }}$ which is represented by some coefficients $A^{\mathbf {\alpha }}$:

$$\begin{aligned} \varPsi = \sum _{\mathbf {\alpha } \in I} A^{\mathbf {\alpha }} \varphi ( {\theta }^{\alpha }, \cdot ), \end{aligned}$$

yields new coefficients $B^{\mathbf {\alpha }}$:

$$\begin{aligned} \begin{pmatrix} B^{(0,0)} \\ B^{(1,1)} \\ B^{(1,0)} \\ B^{(0,1)} \end{pmatrix} = \hat{T}_h^{2h} ( {\theta },\sigma ) \begin{pmatrix} A^{(0,0)} \\ A^{(1,1)} \\ A^{(1,0)} \\ A^{(0,1)} \end{pmatrix}. \end{aligned}$$

$\hat{T}_h^{2\,h} ( {\theta },\sigma ) \in \mathbb {C}^{4 \times 4}$ is a matrix depending on the frequency ${\theta }$ and the exponent $\sigma$. It represents the operator $T_h^{2h}$ with respect to $E_h^{{\theta }}$. In order to determine the convergence behavior of the twogrid method, we study how the coefficients $A^{\mathbf {\alpha }}$ are transformed by $\hat{T}_h^{2h} ( {\theta },\sigma )$ into new coordinates $B^{\mathbf {\alpha }}$. Therefore, the local convergence factor $\rho _{loc}( \sigma )$ depending on the complex shift exponent in the twogrid operator is taken into account:

$$\begin{aligned} \rho _{loc}( \sigma ) = \sup \left\{ \left. \rho ( \hat{T}_h^{2h}( {\theta },\sigma ) ) \right| {\theta } \in T^{low}, {\theta } \notin \varLambda \right\} . \end{aligned}$$

Thereby $\varLambda \subset \left( -\pi ,\pi \right] ^2$ denotes the set of frequencies for which $\hat{L}_{2\,h}$ and $\hat{L}_{h}$ are not invertible and $\rho ( \hat{T}_h^{2h}( {\theta },\sigma ) )$ is the spectral radius of $\hat{T}_h^{2h}( {\theta },\sigma )$. This means that finding the convergence factor for a fixed frequency is reduced to finding the spectral radius of a $4\times 4$ matrix. To obtain the convergence rate of the twogrid solver, one has to search for the maximum of this quantity in the low frequency space $T^{low}$. The search for the maximal spectral radius is restricted to $T^{low}$, since the pattern of $\rho ( \hat{T}_h^{2h} )$ in $T^{low}$ is extended periodically to the whole frequency domain $\varTheta$. Furthermore, an important feature of a multigrid solver is the damping of the low frequency error components [41].

To compute $\hat{T}_h^{2h} ( {\theta },\sigma )$, the operators occurring in the definition of $T_h^{2h}$ are represented with respect to the harmonic space $E_h^{{\theta }}$ by means of $4\times 4$ matrices. Exceptions are the prolongation, restriction and $L^{-1}_{2h}$ operator whose matrices are given by $1\times 4$, $4\times 1$ and $1\times 1$ matrices, since these are mappings related to $E_h^{{\theta }}$ and $E_{2h}^{{\theta }}$. The latter space is the space of harmonics with respect to the coarse grid $\textbf{G}_{2h}$:

$$\begin{aligned} E_{2h}^{{\theta }} = \text {span}\left\{ \left. \varphi _{2h}( {\theta }^{\alpha }, \cdot ) \right| \alpha \in \left\{ (0,0) \right\} \right\} . \end{aligned}$$

The matrices defined with respect to the harmonic spaces are indicated by a hat:

$$\begin{aligned} \hat{T}_h^{2h}( {\theta },\sigma )&= \hat{S}_h^{\nu _2}( {\theta },\sigma ) \hat{K}_h^{2h}( {\theta },\sigma ) \hat{S}_h^{\nu _1}( {\theta },\sigma ), \\ \ \\ \hat{K}_h^{2h}( {\theta },\sigma )&= \hat{I}_h - \hat{I}_{2h}^h( {\theta } ) \hat{L}_{2h}^{-1}( 2{\theta },\sigma ) \hat{I}_h^{2h}( {\theta } ) \hat{L}_{h}( {\theta },\sigma ). \end{aligned}$$

In the following, we list the matrices with respect to the harmonic spaces. Obviously, the matrix for the identity operator $\hat{I}_h$ is given by the identity matrix. The discretization operator $L_h$ can also be represented by a diagonal matrix [12] [41, Chapter 2]:

$$\begin{aligned} \hat{L}_h( {\theta },\sigma ) = \text {diag} \begin{pmatrix} \tilde{L}_h( {\theta }^{(0,0)}, \sigma ),&\tilde{L}_h( {\theta }^{(1,1)}, \sigma ),&\tilde{L}_h( {\theta }^{(1,0)}, \sigma ),&\tilde{L}_h( {\theta }^{(0,1)}, \sigma ) \end{pmatrix}. \end{aligned}$$

The symbol of $L_h$ is given by:

$$\begin{aligned} \tilde{L}_h( {\theta }, \sigma )&= \frac{1}{3} ( 8- 2 \cos (\theta _1) - 2 \cos (\theta _2) - 4 \cos (\theta _1) \cos (\theta _2) ) \\&\quad -\frac{\lambda h^2}{36} (16+8\cos (\theta _1) +8\cos (\theta _2) + 4 \cos (\theta _1) \cos (\theta _2) ). \end{aligned}$$

The matrix for the coarse grid operator and ${\theta } = {\theta }^{(0,0)}$ is given by:

$$\begin{aligned} \hat{L}_{2h}({\theta },\sigma ) = \tilde{L}_{2h}({\theta },\sigma ),\; L_{2h}(\sigma )\varphi _{2h}( 2{\theta },\textbf{x} ) = \tilde{L}_{2h}({\theta },\sigma ) \varphi _{2h}( 2{\theta },\textbf{x} ). \end{aligned}$$

The exact formula for $\tilde{L}_{2h}({\theta },\sigma )$ reads as follows:

$$\begin{aligned} \tilde{L}_{2h}(2{\theta },\sigma )&= \frac{1}{3} ( 8- 2 \cos (2\theta _1) - 2 \cos (2\theta _2) - 4 \cos (2\theta _1) \cos (2\theta _2) ) \\&\quad -\frac{\lambda h^2}{9} (16+8\cos (2\theta _1) +8\cos (2\theta _2) + 4 \cos (2\theta _1) \cos (2\theta _2) ). \end{aligned}$$

In a next step the matrices for the restriction operator and the prolongation operator are derived:

$$\begin{aligned} \hat{I}_{2h}^h ({\theta } ) = \begin{pmatrix} \tilde{I}_{2h}^h( {\theta }^{(0,0)} ) \\ \tilde{I}_{2h}^h( {\theta }^{(1,1)} ) \\ \tilde{I}_{2h}^h( {\theta }^{(1,0)} ) \\ \tilde{I}_{2h}^h( {\theta }^{(0,1)} ) \end{pmatrix} = \begin{pmatrix} (1+\cos (\theta _1) ) (1+\cos (\theta _2) ) \\ (1-\cos (\theta _1) ) (1-\cos (\theta _2) ) \\ (1-\cos (\theta _1) ) (1+\cos (\theta _2) ) \\ (1+\cos (\theta _1) ) (1-\cos (\theta _2) ) \end{pmatrix} \text { and } \hat{I}_h^{2h}({\theta } ) = \frac{1}{4} \cdot (\hat{I}_{2h}^h ({\theta } ) )^T. \end{aligned}$$

Thereby the matrix entries are determined by the following equation:

$$\begin{aligned} I_{2h}^h \varphi _{2h}( 2{\theta },\textbf{x} ) = \sum _{\alpha } \tilde{I}_{2h}^h( {\theta }^\alpha ) \varphi _{h}( {\theta }^\alpha ,\textbf{x} ) \text { and } I_h^{2h} \varphi _{h}( {\theta }^\alpha ,\textbf{x} ) = \tilde{I}_h^{2h}({\theta }^\alpha ) \varphi _{2h}( 2{\theta }^{(0,0)},\textbf{x} ),\; \alpha \in I. \end{aligned}$$

Fig. 11

Convergence factors $\rho$ in the frequency space $T^{low}$ in case of four different shifts $\sigma \in \left\{ 1.0,1.25,1.75,2.0 \right\}$

Fig. 12

Convergence factor $\rho$ for $\theta _2=0$ and $\theta _1 \in \left[ -\frac{\pi }{2},\frac{\pi }{2} \right]$.The exponents for the shift are chosen as follows: $\sigma \in \left\{ 1.00,1.25,1.50,1.75,2.00 \right\}$

We point out that the restriction and prolongation operator are independent of k and $\sigma$. Summarizing the previous considerations, one obtains the matrix $\hat{K}_h^{2h}({\theta },\sigma )$. Finally, it remains to specify the matrix for the smoother. According to [41, Chapter 4], it is given by a diagonal matrix:

$$\begin{aligned} \hat{S}_{h}( {\theta },\sigma ) = \text {diag} \begin{pmatrix} \tilde{S}_{h}( {\theta }^{(0,0)},\sigma ),&\tilde{S}_{h}( {\theta }^{(1,1)},\sigma ),&\tilde{S}_{h}( {\theta }^{(1,0)},\sigma ),&\tilde{S}_{h}( {\theta }^{(0,1)},\sigma ) \end{pmatrix}. \end{aligned}$$

A straightforward calculation yields for ${\theta } = {\theta }^{(0,0)}$:

$$\begin{aligned} \tilde{S}_{h}( {\theta }, \sigma )&= (1-\omega ) + \omega \frac{\frac{2}{3} + \lambda h^2 \frac{8}{36}}{\frac{8}{3} - \lambda h^2 \frac{16}{36}} \cos ( \theta _1 ) \\&\quad + \omega \frac{\frac{2}{3} + \lambda h^2 \frac{8}{36}}{\frac{8}{3} - \lambda h^2 \frac{16}{36}} \cos ( \theta _2 ) + \omega \frac{\frac{4}{3} + \lambda h^2 \frac{4}{36}}{\frac{8}{3} - \lambda h^2 \frac{16}{36}} \cos ( \theta _1 ) \cos ( \theta _2 ). \end{aligned}$$

Finally, we have all ingredients at hand to assemble the matrix $\hat{T}_h^{2h}$ so that $\rho _{loc}( \sigma )$ can be computed numerically.

In Fig. 11, we show some typical surfaces for $\rho$ in the space consisting of low frequencies for $\nu _1=\nu _2=3$, $\omega = 2/3$, $h=2^{-5}$ and $k=0.35/h$. It can be observed that the maximum of the spectral radius $\rho ( \hat{T}_h^{2h} )$, i.e., $\rho _{loc}$ increases as $\sigma$ is decreased from $\sigma = 2$ towards $\sigma = 1$. In addition to that, there is a symmetry with respect to the axes of the coordinate system. This can be explained by the fact that the entries of the Fourier symbols are composed of cosine functions (see “Appendix A”). We further observe numerically that the maxima are located on the axes for $\theta _1 = 0$ or $\theta _2 =0$ (see Fig. 1). Motivated by these observations, we restrict our search for the maximal spectral radius of $\hat{T}_h^{2h}$ to the line with $\theta _2 =0$ (see Fig. 12).