Parallel Controllability Methods For the Helmholtz Equation

TL;DR

This paper introduces parallel controllability methods for solving high-frequency Helmholtz equations by transforming the problem into the time domain, resulting in scalable algorithms that handle large-scale problems efficiently.

Contribution

The paper develops robust, parallel controllability algorithms for the Helmholtz equation using first and second-order wave formulations, applicable to general boundary-value problems.

Findings

Achieves high accuracy and convergence in Helmholtz solutions

Demonstrates strong scalability on massively parallel architectures

Handles problems with up to a billion unknowns

Abstract

The Helmholtz equation is notoriously difficult to solve with standard numerical methods, increasingly so, in fact, at higher frequencies. Controllability methods instead transform the problem back to the time-domain, where they seek the time-harmonic solution of the corresponding time-dependent wave equation. Two different approaches are considered here based either on the first or second-order formulation of the wave equation. Both are extended to general boundary-value problems governed by the Helmholtz equation and lead to robust and inherently parallel algorithms. Numerical results illustrate the accuracy, convergence and strong scalability of controllability methods for the solution of high frequency Helmholtz equations with up to a billion unknowns on massively parallel architectures.

Tables

Table 1. Table 1: 2D-Marmousi model: P 2 superscript 𝑃 2 P^{2} -FE with 15 15 15 points per wave length

Frequency	Wave number	#Unknowns	#Nodes
$\nu$ [Hz]	$k=\omega /c=2\pi \nu /c$	$ndof$	24 cores per node
$10$	$11$ – $42$	$1^{\prime }{658}^{\prime }443$	$1$–$8$
$20$	$22$ – $84$	$6^{\prime }{628}^{\prime }881$	$1$–$16$
$40$	$45$ – $168$	${26}^{\prime }{505}^{\prime }761$	$8$–$64$
$60$	$68$ – $252$	${59}^{\prime }{630}^{\prime }641$	$16$–$128$
$80$	$91$ – $336$	${106}^{\prime }{003}^{\prime }521$	$16$–$128$
$160$	$182$ – $671$	${423}^{\prime }{975}^{\prime }041$	$64$–$256$
$250$	$285$ – $1048$	$1^{\prime }{035}^{\prime }{241}^{\prime }009$	$128$–$512$

Table 2. Table 2: 3D-cavity: CMCG methods with P 1 superscript 𝑃 1 P^{1} -FEM. As η 𝜂 \eta increases, the ratio h ​ k 3 / 2 ℎ superscript 𝑘 3 2 hk^{3/2} remains constant to avoid pollution errors [ 30 ] .

Frequency	#Unknowns	#Tetrahedra	CG iterations	#Nodes
$\nu =2\pi \omega$	$ndof$			24 cores per node
$2$	$8.17\cdot {10}^5$	$5^{\prime }{051}^{\prime }049$	$239$	$1$–$8$
$3$	$5.22\cdot {10}^6$	${31}^{\prime }{190}^{\prime }000$	$440$	$2$–$32$
$4$	$1.9\cdot {10}^7$	${114}^{\prime }{391}^{\prime }112$	$607$	$32$–$96$
$6$	$1.18\cdot {10}^8$	${703}^{\prime }{590}^{\prime }464$	$578$	$64$–$128$

Equations

$$-\Delta u(x)-k^2(x)u(x)$$

$$\frac{\partial u(x)}{\partial n}-ik(x)u(x)$$

$$\frac{\partial u(x)}{\partial n}$$

$$u(x)$$

$$\frac{1}{c^2(x)}\frac{{\partial }^{2}y(x,t)}{{\partial }^{2}t}-\Delta y(x,t)$$

$$\frac{\partial y(x,t)}{\partial n}+\frac{1}{c(x)}\frac{\partial y(x,t)}{\partial t}$$

$$\frac{\partial y(x,t)}{\partial n}$$

$$y(x,t)$$

$$y(x,0)=v_0(x),\frac{\partial y(x,0)}{\partial t}$$

$$u=v_0+\frac{i}{\omega }{v}_1,v_0,v_1\in H^1(\Omega ).$$

$$J(v_0,v_1)=\frac{1}{2}{\int }_{\Omega }|\nabla y(x,T)-\nabla v_0(x){|}^{2}dx+\frac{1}{2}{\int }_{\Omega }\frac{1}{c^2(x)}(y_t(x,T)-v_1(x){)}^{2}dx,$$

$$(y(\cdot ,t),\phi )=(u{e}^{-i\omega t},\phi )+(\lambda +\eta t,\phi )+\sum _{|\ell |>1}({\gamma }_{\ell }{e}^{i\omega \ell t},\phi )$$

$$-\Delta {\gamma }_{\ell }(x)-(\ell k(x){)}^2{\gamma }_{\ell }(x)$$

$$\frac{\partial {\gamma }_{\ell }(x)}{\partial n}+i\ell k(x){\gamma }_{\ell }(x)$$

$$\frac{\partial {\gamma }_{\ell }(x)}{\partial n}$$

$${\gamma }_{\ell }(x)$$

$$(v,\phi )=(u,\phi )+(\lambda +\frac{i}{\omega }\eta ,\phi )+\sum _{|\ell |>1}({\alpha }_{\ell }+i\ell {\beta }_{\ell },\phi ),\forall \phi \in H_D^1.$$

$$\lambda=\frac{1}{\|k\|_{L^{2}(\Omega)}^{2}+i|k|_{L^{1}(\Gamma_{S})}}\bigg{(}\int_{\Omega}k^{2}v+i\int_{\Gamma_{S}}kv+\int_{\Omega}f+\int_{\Gamma_{S}}g_{S}+\int_{\Gamma_{N}}g_{N}\bigg{)}.$$

$$\widehat{y}(x):=\frac{1}{T}\int_{0}^{T}\big{(}y(x,t)+\frac{i}{\omega}y_{t}(x,t)\big{)}\operatorname{e}^{i\omega t}\,dt.$$

$$\hat{y}(x)$$

$$=\frac{1}{T}{\int }_0^{T}u{e}^{-i\omega t}{e}^{i\omega t}dt-\frac{i\eta }{\omega }=u-\frac{i\eta }{\omega }.$$

$$u(x)=\hat{y}(x)+\frac{i\eta }{\omega },x\in \Omega$$

$$-{\int }_{\Omega }{k}^2(x)u(x)dx={\int }_{\Omega }f(x)dx+{\int }_{\partial \Omega }{g}_N(x)ds.$$

$$\frac{i\eta }{\omega }$$

$$⟨J^{\prime }(v),\delta v⟩$$

$$\frac{1}{c^2(x)}\frac{{\partial }^2}{{\partial }^{2}t}p(x,t)-\Delta p(x,t)$$

$$\frac{\partial p(x,t)}{\partial n}-\frac{1}{c}\frac{\partial }{\partial t}p(x,t)$$

$$\frac{\partial p(x,t)}{\partial n}$$

$$p(x,t)$$

$$p(x,T)=p_0(x),\frac{\partial p(x,T)}{\partial t}$$

$$p_0(x)$$

$${\int }_{\Omega }\frac{p_1(x)}{c^2(x)}w(x)dx$$

$$(\nabla {\tilde{g}}_0,\nabla \phi )$$

$$={\int }_{\Omega }\nabla (v_0(x)-y(x,T))\cdot \nabla \phi (x)-\frac{1}{c^2(x)}{p}_t(x,0)\phi (x)dx$$

$$+{\int }_{{\Gamma }_S}\frac{1}{c(x)}p(x,0)\phi (x)ds,\forall \phi \in H_D^1,$$

$${\tilde{g}}_1$$

$$\sqrt{\frac{\Vert \nabla r_0^{(\ell +1)}{\Vert }_{L^2(\Omega )}^2+\Vert (1/c)r_1^{(\ell +1)}{\Vert }_{L^2(\Omega )}^2}{\Vert \nabla r_0^{(0)}{\Vert }_{L^2(\Omega )}^2+\Vert (1/c)r_1^{(0)}{\Vert }_{L^2(\Omega )}^2}}\le tol.$$

$$u_h=v_0^{(\ell )}+\frac{i}{\omega }{v}_1^{(\ell )}.$$

$$\frac{1}{c^2(x)}{v}_t(x,t)-\nabla \cdot \mathbf{p}(x,t)$$

$$\frac{\partial }{\partial t}\mathbf{p}(x,t)$$

$$\mathbf{p}(x,t)\cdot \mathbf{n}+\frac{1}{c(x)}v(x,t)$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

The efficient numerical solution of the Helmholtz equation is fundamental to the simulation of time-harmonic wave phenomena in acoustics, electromagnetics or elasticity. As the time frequency $\omega>0$ increases, so does the size of the linear system resulting from any numerical discretization in order to resolve the increasingly smaller wave lengths. With the increase in frequency, however, the performance of standard preconditioners based on multigrid, incomplete factorization or domain decomposition approaches, originally developed for positive definite Laplace-like equations, rapidly deteriorates [1].

In recent years, a growing number of increasingly sophisticated preconditioners has been proposed for the iterative solution of the Helmholtz equation; ”Shifted Laplacian” preconditioners [2], for instance, have led to modern multigrid [3, 4] and domain decomposition preconditioners [5, 6]. While some of those preconditioners may achieve a desirable frequency independent convergence behavior in special situations [7], that optimal behavior is often lost in the presence of strong heterogeneity. Moreover, they are typically tied to a special discretization or fail to achieve optimal scaling on parallel architectures.

Controllability methods (CM) offer an alternative approach for the numerical solution of the Helmholtz equation. Instead of solving the problem directly in the frequency domain, we first transform it back to the time domain where we seek the corresponding time-dependent periodic solution, $y(\cdot,t)$, with known period $T=2\pi/\omega$. By minimizing an energy functional $J(v_{0},v_{1})$ which penalizes the mismatch after one period, controllability methods iteratively adjust the (unknown) initial condition $(v_{0},v_{1})$ thereby steering $y(\cdot,t)$ towards the desired periodic solution. Once the minimizer of $J$ has been found, we immediately recover from it the solution of the Helmholtz equation. As the CM combines the numerical integration of the time-dependent wave equation with a conjugate gradient (CG) iteration, it is remarkably robust and inherently parallel.

In [8], Bristeau et al. proposed the first CM for sound-soft scattering problems based on the wave equation in standard second-order form. Since the initial condition $(v_{0},v_{1})$ then lies in $H^{1}\times L^{2}$, the original formulation requires the solution of a coercive elliptic problem at each CG iteration. Heikkola et al. in [9, 10] presented a higher-order version by using spectral FE and the classical fourth-order Runge-Kutta (RK) method. For more general boundary-value problems, such as wave scattering from sound-hard obstacles, inclusions, or wave propagation in physically bounded domains, the original CM will generally fail because the minimizer of $J$ is no longer unique. In [11], we proposed alternative energy functionals which restore uniqueness, albeit at a small extra computational cost, for general boundary-value problems governed by the Helmholtz equation.

More recently, Glowinski and Rossi [12] proposed a CM based on the wave equation in first-order (or mixed) form using classical Raviart-Thomas (RT) finite elements. As $(v_{0},v_{1})$ then lies in $L^{2}\times(L^{2})^{d}$, the solution of an elliptic problem at each CG iteration is no longer necessary and the CM becomes in principle trivially parallel. Still, the lack of availability of mass-lumping for RT elements again nullifies the main advantage of the first-order formulation because the mass-matrix now needs to be ”inverted” at each time-step.

Here we revisit the original CM from [8, 12] and consider two distinct discretizations, which both lead to highly efficient and inherently parallel methods. In Section 2, we recall the CMCG method based on the wave equation in second-order form and propose a filtering procedure which permits the use of the original energy functional $J$, regardless of the boundary conditions. Next, in Section 3, we consider the CM based on the wave equation in first-order form and again show how to extend it to arbitrary boundary-value problems governed by the Helmholtz equation. Thanks to a recent hybrid discontinuous Galerkin (HDG) method [13], which automatically yields a block-diagonal mass-matrix, the time integration of the wave equation then becomes truly explicit and the entire CMCG approach trivially parallel. In Section 4, we perform a series of numerical experiments to illustrate the accuracy, convergence behavior and inherent parallelism of the CMCG approach. In particular, we apply it to large-scale high-frequency Helmholtz problems with up to a billion unknowns to demonstrate its strong scalability on massively parallel architectures.

2 Controllability methods for the second-order formulation

2.1 Time-harmonic waves

We consider a time-harmonic wave field $u(x)$ in a bounded connected computational domain $\Omega\subset\mathbb{R}^{d}$, $d\leq 3$, with a Lipschitz boundary $\Gamma$. The boundary consists of three disjoint components, $\Gamma=\Gamma_{D}\cup\Gamma_{N}\cup\Gamma_{S}$ where we impose a Dirichlet, Neumann and impedance (or Sommerfeld-like absorbing) boundary condition, respectively; the boundary condition is omitted whenever the corresponding component is empty. In $\Omega$, the wave field $u$ hence satisfies the Helmholtz equation

[TABLE]

where $\omega>0$ is the (angular) frequency, $c(x)>0$ the wave speed, $k(x)=\omega/c(x)$ the wave number, $n$ the unit outward normal, and $f$, $g_{N}$, $g_{S}$ and $g_{D}$ are known and may vanish.

The above formulation is rather general and encompasses most common applications such as sound-soft scattering problems with $\Gamma_{S}\neq\emptyset$ and $\Gamma_{D}\neq\emptyset$, sound-hard scattering problems with $\Gamma_{S}\neq\emptyset$ and $\Gamma_{N}\neq\emptyset$, or physically bounded domains with $\Gamma_{S}=\emptyset$. We shall always assume for any particular choice of $\omega$, $c(x)$, or combination of boundary conditions that (2.1) is well-posed and has a unique solution $u\in H^{1}(\Omega)$.

Instead of solving the Helmholtz equation directly in the frequency domain, we now reformulate (2.1) in the time domain. Then, the corresponding time-harmonic wave field, $u(x)\operatorname{e}^{-i\omega t}$, satisfies the (real-valued) time-dependent wave equation

[TABLE]

for the (unknown) initial values $v_{0}=u$ and $v_{1}=\omega\operatorname{Im}\left\{u\right\}$.

For sound-soft scattering problems (2.1), where $|\Gamma_{D}|>0$ and $|\Gamma_{S}|>0$, Bristeau et al. [8, 14] proposed to determine $u(x)$ via controllability by computing a time-periodic solution $y(x,t)$ of (2.2) with period $T={2\pi}/{\omega}$. Once the initial values $v_{0},v_{1}$ of $y$ are known, the solution $u$ of the original Helmholtz equation (2.1) is immediately given by

[TABLE]

To determine $v_{0}$ and $v_{1}$, the problem is reformulated as a least-squares optimization problem over $H^{1}(\Omega)\times L^{2}(\Omega)$ for the quadratic cost functional

[TABLE]

where $y$ satisfies (2.2) with the initial values $v_{0}$ and $v_{1}$. The functional $J$ measures in the energy norm the mismatch between the solution of (2.2) at the initial time and after one period. It is non-negative and convex, while $J(v_{0},v_{1})=0$ if, and only if, $\nabla y(\cdot,T)=\nabla y(\cdot,0)$ and $y_{t}(\cdot,T)=y_{t}(\cdot,0)$ for any given initial values $(v_{0},v_{1})$; in particular, $J(v_{0},v_{1})=0$ if $v_{0}=u$ and $v_{1}=\omega\operatorname{Im}\left\{u\right\}$.

For more general scattering problems, however, $J$ is no longer strictly convex as the $T$-periodicity of $y_{t}$ and $\nabla y$ no longer guarantees a unique periodic solution $y$ of (2.2). Instead, for the general boundary-value problem (2.1), the situation is more complicated and summarized in the following theorem [11] .

Theorem 1**.**

Let $u\in H^{1}(\Omega)$ be the unique solution of (2.1) and $y\in C^{0}([0,T];H^{1}(\Omega))\cap C^{1}([0,T];L^{2}(\Omega))$ be a (real-valued) solution of (2.2) with initial values $(v_{0},v_{1})\in H^{1}(\Omega)\times L^{2}(\Omega)$. If $\nabla y$ and $y_{t}$ are time periodic with period $T=2\pi/\omega$, then $y$ admits the Fourier series representation

[TABLE]

for any $\varphi\in H_{D}^{1}$, where the constants $\lambda,\eta\in\mathbb{R}$ and the eigenfunctions $\gamma_{\ell}={\alpha}_{\ell}+i{\beta}_{\ell}$, ${\alpha}_{\ell},{\beta}_{\ell}\in H^{1}(\Omega)$, $|\ell|>1$ satisfy

[TABLE]

Let $v=v_{0}+({i}/{\omega})\ v_{1}$. Then $v$ satisfies

[TABLE]

*Furthermore, if $|\Gamma_{S}|>0$, then $\eta=0$. If $|\Gamma_{D}|>0$, then $\lambda=\eta=0$.

Here H_{D}^{1}:=\{w\in H^{1}(\Omega):w=0\text{ on \Gamma_{D}}\} and $(\cdot,\cdot)$ denotes the standard $L^{2}(\Omega)$ inner product.*

Proof.

See [11]. ∎

For sound-soft scattering problems ($|\Gamma_{S}|>0,|\Gamma_{D}|>0$), where both Dirichlet and Sommerfeld-like absorbing boundary conditions are imposed on $\Gamma$, all the eigenfunctions $\gamma_{\ell}$, $|\ell|>1$, and the constants $\lambda$, $\eta$ in (2.7) vanish identically. Thus, the minimizer $v=v_{0}+({i}/{\omega})v_{1}$ of $J$ in (2.4) then coincides with $u$.

For scattering problems from sound-hard obstacles or penetrable inclusions ($|\Gamma_{S}|>0$, $|\Gamma_{D}|=0$), the eigenfunctions $\gamma_{\ell}$ and the constant $\eta$ in (2.7) still vanish identically, yet the constant $\lambda$ may be nonzero. Given any minimizer $v=u+\lambda$ of $J$, we can recover $u$ by subtracting the spurious shift $\lambda$ using the compatibility condition:

[TABLE]

In fact, any impedance condition (2.1b) that includes a positive (or negative) definite zeroth order term, such as a more accurate absorbing boundary condition [15, 16], also circumvents the indeterminacy due to $\lambda$. For wave propagation in physically bounded domains ($|\Gamma_{S}|=0$), the eigenfunctions $\gamma_{\ell}$ and the constants $\lambda,\eta$ in (2.7) typically do not vanish. However, we can always restore uniqueness by replacing $J$ with an alternative energy functional, thereby incurring a small increase in computational cost – see [11].

2.2 Fundamental frequency extraction via filtering

From Theorem 1 we conclude that a minimizer of $J$ generally yields a time-dependent solution $y$ of (2.2), which contains a constant shift determined by $\lambda$, a linearly growing part determined by $\eta$, and higher frequency harmonics determined by $\gamma_{\ell}$, all superimposed on the desired time-harmonic field $u$ with fundamental frequency $\omega$. Those spurious modes can be eliminated by replacing $J$ with an alternative energy functional at a small extra computational cost [11]. Instead we now propose an alternative approach via filtering which removes all spurious modes without requiring a modified energy functional.

Let $y(x,t)$ be the time-dependent solution of (2.2) that corresponds to a minimizer $(v_{0},v_{1})$ of $J$. Next, we define \widehat{y}\in\{w\in H^{1}(\Omega)\ |\ w=g_{D}\text{ on \Gamma_{D}}\} as

[TABLE]

To extract $u(x)$ from $y(x,t)$, we now take advantage of the mutual orthogonality of different time harmonics $\exp(i\omega\ell t)$ in $L^{2}(0,T)$. Hence, we multiply (2.5) with $\operatorname{e}^{i\omega t}$ and integrate in time over $(0,T)$ to obtain

[TABLE]

This yields

[TABLE]

where $\lambda$ and all $\gamma_{\ell}$ have vanished but the constant $\eta$ is still undetermined.

If $|\Gamma_{S}|>0$ or $|\Gamma_{D}|>0$, Theorem 1 implies that $\eta=0$ and thus $u(x)=\widehat{y}(x)$. Otherwise in the pure Neumann case ($\Gamma=\Gamma_{N}$), we determine $\eta$ by integrating (2.10), multiplied by $k^{2}(x)$, over $\Omega$ and using the compatibility condition

[TABLE]

from (2.1a). This immediately yields the remaining constant

[TABLE]

We summarize the above derivation in the following proposition.

Proposition 1**.**

Let $u\in H^{1}(\Omega)$ be the unique solution of (2.1) and $y$ the time dependent solution of (2.2) corresponding to a minimizer $(v_{0},v_{1})\in H^{1}(\Omega)\times L^{2}(\Omega)$ of $J$, i.e. $J(v_{0},v_{1})=0$. Then $u$ is given by (2.10) with $\eta=0$ if $|\Gamma_{S}|>0$ or $|\Gamma_{D}|>0$, and with $\eta$ given by (2.12) when $\Gamma_{N}=\partial\Omega$.

Not only does the above filtering approach allow us to use the original cost functional $J$, it also involves a negligible computational effort or storage amount, as the time integral for $\widehat{y}$ can be calculated cumulatively via numerical quadrature during the solution of the wave equation (2.2).

2.3 The CMCG Algorithm

To minimize the quadratic cost functional $J$ defined by (2.4) over $H^{1}(\Omega)\times L^{2}(\Omega)$, a natural choice is the conjugate gradient (CG) method [8], which requires the Fréchet derivative of $J$ at $v=(v_{0},v_{1})$:

[TABLE]

Here $\delta v=(\delta v_{0},\delta v_{1})$ denotes an arbitrary perturbation, $\langle\cdot,\cdot\rangle$ the standard duality pairing, and $p$ the solution of the adjoint (backward) wave equation:

[TABLE]

The derivation of (2.3) and (2.15) can be found in [8]. In each CG iteration the derivative $J^{\prime}(v)$ requires the solution of the forward and backward (adjoint) wave equations (2.2) and (2.15) over one period $[0,T]$. Moreover, each CG iteration requires an explicit (Riesz) representer $\tilde{g}=(\tilde{g}_{0},\tilde{g}_{1})\in H_{D}^{1}(\Omega)\times L^{2}(\Omega)$ of the gradient $g=(g_{0},g_{1})=J^{\prime}(v)$ defined in (2.3), which is determined by solving the symmetric and coercive elliptic problem [8, 17]:

[TABLE]

For the sake of completeness, we list the full CMCG Algorithm – see [8, 11]:

**CMCG Algorithm. ** **

1. (1)

Initialize $v^{(0)}=(v^{(0)}_{0},v^{(0)}_{1})$ (initial guess). 2. (2)

Solve the forward and the backward wave equations (2.2) and (2.15) to determine the gradient of $J$, $g^{(0)}=J^{\prime}(v^{(0)})$, defined by (2.3). 3. (3)

Solve the coercive elliptic problem (2.16) with $g=g^{(0)}$ to determine the new search direction $\tilde{g}^{(0)}$. 4. (4)

Set $r^{(0)}=d^{(0)}=\tilde{g}^{(0)}$. 5. (5)

For $\ell=1,2,\ldots$

1. 5.1 

Solve the wave equation (2.2) with $f=g_{D}=g_{S}=g_{N}=0$ and the initial values $d^{(\ell)}=(d^{(\ell)}_{0},d^{(\ell)}_{1})$ and the backward wave equation (2.15). Compute the gradient $g^{(\ell)}=J^{\prime}(d^{(\ell)})$ defined by (2.3). 2. 5.2 

Solve the coercive elliptic problem (2.16) with $g=g^{(\ell)}$ to get $\tilde{g}^{(\ell)}$. 3. 5.3 

$\alpha_{\ell}=\dfrac{\|\ \nabla r^{(\ell)}_{0}\|_{L^{2}(\Omega)}^{2}+\|(1/c)\ r^{(\ell)}_{1}\|_{L^{2}(\Omega)}^{2}}{(\ \nabla\tilde{g}^{(\ell)}_{0},\nabla d^{(\ell)}_{0})_{L^{2}(\Omega)}+((1/c^{2})\ \tilde{g}^{(\ell)}_{1},d^{(\ell)}_{1})_{L^{2}(\Omega)}}$** 4. 5.4 

$v^{(\ell+1)}=v^{(\ell)}-\alpha_{\ell}d^{(\ell)}$** 5. 5.5 

$r^{(\ell+1)}=r^{(\ell)}-\alpha_{\ell}\tilde{g}^{(\ell)}$** 6. 5.6 

$\beta_{\ell}=\dfrac{\|\nabla r^{(\ell+1)}_{0}\|_{L^{2}(\Omega)}^{2}+\|(1/c)\ r^{(\ell+1)}_{1}\|_{L^{2}(\Omega)}^{2}}{\|\nabla r^{(\ell)}_{0}\|_{L^{2}(\Omega)}^{2}+\|(1/c)\ r^{(\ell)}_{1}\|_{L^{2}(\Omega)}^{2}}$** 7. 5.7 

$d^{(\ell+1)}=r^{(\ell+1)}+\beta_{\ell}d^{(\ell)}$** 8. 5.8 

Stop when the relative residual lies below the given tolerance $tol$

[TABLE] 6. (6)

Return approximate solution $u_{h}$ of (2.1) given by

[TABLE]

Since $\tilde{g}_{0}\in H^{1}(\Omega)$, the updates of $r_{0}^{(k)}$, $d_{0}^{(k)}$ and $v_{0}^{(k)}$ in Steps 5.4, 5.5 and 5.7 in the CMCG Algorithm also remain in $H^{1}(\Omega)$. We emphasize that (2.16a) is independent of $\omega$ and leads to a symmetric and positive definite linear system, which can be solved efficiently and in parallel with standard numerical (multigrid, domain decomposition, etc.) methods [18, 6].

3 Controllability methods for first-order formulations

The CMCG Algorithm from Section 2.3 iterates on the initial value $(v_{0},v_{1})\in H^{1}(\Omega)\times L^{2}(\Omega)$ of the second-order wave equation (2.2) until its solution is $T$-time periodic. However, the gradient of the cost functional $J(v_{0},v_{1})$, which is needed during the CG update, only lies in the dual space $H^{-1}(\Omega)\times L^{2}(\Omega)$. To ensure that the solution remains sufficiently regular and in $H^{1}(\Omega)\times L^{2}(\Omega)$, the corresponding Riesz representative is computed at every CG iteration by solving the strongly elliptic problem (2.16a). In [12], Glowinski et al. derived an equivalent first-order formulation for sound-soft scattering problems, where the solution instead lies in $(L^{2}(\Omega))^{d+1}$, which is reflexive. As a consequence, all CG iterates automatically lie in the correct solution space $(L^{2}(\Omega))^{d+1}$, while the solution of (2.16a) is no longer needed.

3.1 First-order formulation for general boundary conditions

Again, we always assume for any particular choice of $\omega$, $c(x)$, $f$ and combination of boundary conditions that (2.1) has a unique solution $u\in H^{1}(\Omega)$. Following [12], we now let $v=y_{t}$, ${\bf p}=\nabla y$ and rewrite the time-dependent wave equation (2.2) in first-order form:

[TABLE]

Hence, the solution $({\bf p},v)$ of (3.1) lies in the function space $\mathcal{Q}$ [19, 20],

[TABLE]

In terms of ${\bf p}$ and $v$, the energy functional $J$ defined in (2.4) now becomes

[TABLE]

where $({\bf p},v)$ solves (3.1) with initial value $({\bf p}_{0},v_{0})\in H(\operatorname{div};\Omega)\times L^{2}(\Omega)$.

The CMCG Algorithm for the first-order formulation is identical to that for the second-order formulation from Section 2.3 except for Steps 2 and 5.1, where $J^{\prime}$ is now replaced by $\widehat{J}^{\prime}$:

[TABLE]

Here $(\delta{\bf p}_{0},\delta v_{0})\in{\bf P}\times L^{2}(\Omega)$ denotes an arbitrary perturbation with

[TABLE]

whereas $({\bf p}^{*},v^{*})\in{\bf P}\times L^{2}(\Omega)$ solves the backward (adjoint) wave equation in first-order form [12], that is (3.1) with $f\equiv g_{S}\equiv g_{N}\equiv g_{D}\equiv 0$ and

[TABLE]

For sound-soft scattering problems ($|\Gamma_{D}|,|\Gamma_{S}|>0$), the functional $\widehat{J}$ always has a unique (global) minimizer, which therefore coincides with the (unique) time-harmonic solution $u(x)\operatorname{e}^{-i\omega t}$ of (3.1). For more general boundary value problems, however, the minimizer of $\widehat{J}$ is not necessarily unique, as shown in the following theorem.

Theorem 2**.**

Let $u\in H^{1}(\Omega)$ be the unique solution of (2.1) and $({\bf p},v)\in\mathcal{Q}$ be a real-valued solution of (3.1) with initial values $({\bf p}_{0},v_{0})\in H(\operatorname{div};\Omega)\times L^{2}(\Omega)$. If ${\bf p}$ and $v$ are time periodic with period $T=2\pi/\omega$, then ${\bf p}$ and $v$ admit the Fourier series representation

[TABLE]

where the constant $\eta\in\mathbb{R}$, $\boldsymbol{\lambda}\in{\bf P}$ with

[TABLE]

and the complex-valued eigenfunctions $\boldsymbol{\gamma}_{\ell}^{p}\in{\bf P}$, $\gamma_{\ell}^{v}\in L^{2}(\Omega)$, $|\ell|>1$ satisfy

[TABLE]

Furthermore, if $|\Gamma_{S}\cup\Gamma_{D}|>0$, then $\eta=0$.

Proof.

Let

[TABLE]

Then $w$ and ${\bf q}$ satisfy (3.1) with $f\equiv g_{D}\equiv g_{S}\equiv g_{N}\equiv 0$ and initial values

[TABLE]

Since ${\bf p}$ and $v$ are $T$-periodic, so are ${\bf q}$ and $w$. Moreover, the mappings

[TABLE]

are $T$-periodic and continuous for any $(\boldsymbol{\psi},\varphi)\in{\bf P}\times L^{2}(\Omega)$ [19]. Hence, they admit the Fourier series representation,

[TABLE]

where $\boldsymbol{\gamma}_{\ell}^{p}\in\mathbb{C}^{d}$ , $\widehat{\gamma}_{\ell}^{v}\in\mathbb{C}$. Next, we define

[TABLE]

which implies that

[TABLE]

We shall now show that $\boldsymbol{\gamma}_{\ell}^{p}$ and $\gamma_{\ell}^{v}$ satisfy (3.8) for all $|\ell|\geq 1$. First, integration by parts, (3.1a)-(3.1b) and the periodicity of ${\bf q}$ and $w$ imply

[TABLE]

Together with definition (3.9) of ${\boldsymbol{\gamma}}_{\ell}^{p}$ and $\gamma_{\ell}^{v}$, we thus immediately obtain

[TABLE]

Since $w(x,t)=0$ for $x\in\Gamma_{D}$, we infer from (3.9) that

[TABLE]

and hence $\gamma_{\ell}^{v}$ satisfies (3.8e). Similarly, (3.8c), (3.8d) follow from the fact that ${\bf q}$ and $w$ satisfy (3.1c), (3.1d) with $g_{N}\equiv g_{S}\equiv 0$. Hence $\boldsymbol{\gamma}_{\ell}^{p}$, $\gamma_{\ell}^{v}$ satisfy (3.8) for all $|\ell|\geq 1$. In fact for $\ell=1$, (3.8) corresponds to (2.1) in first-order formulation with $\boldsymbol{\gamma}_{1}^{p}=\nabla\overline{u}$, $\gamma_{1}^{v}=i\omega\overline{u}$, homogeneous boundary conditions and no sources. By uniqueness, $\boldsymbol{\gamma}_{1}^{p}$ and $\gamma_{1}^{v}$, together with their complex conjugates, are therefore identically zero.

Next, we consider $\boldsymbol{\gamma}_{0}^{p}$, $\gamma_{0}^{v}$. Again, since ${\bf q}$ and $w$ satisfy (3.1a)-(3.1e) with $f=0$ and homogeneous boundary conditions, we obtain from (3.9) with $\ell=0$ and the periodicity of ${\bf q}$ and $w$

[TABLE]

In particular, (3.10)-(3.11) implies with $\varphi=\gamma_{0}^{v}$ and $\boldsymbol{\psi}=\boldsymbol{\gamma}_{0}^{p}$ that

[TABLE]

and hence, $\boldsymbol{\gamma}_{0}^{p}\cdot{\bf n}=0$ on $\Gamma_{S}$, since $c>0$. Moreover, Green’s formula, together with (3.10) and the homogeneous boundary conditions, implies that

[TABLE]

and therefore $\boldsymbol{\lambda}=\boldsymbol{\gamma}_{0}^{p}$ satisfies (3.7).

To show that $\gamma_{0}^{v}$ is constant, we now let $\varphi\in\mathcal{C}_{c}^{\infty}(\Omega)$ and $\boldsymbol{\psi}={\bf e}_{j}\varphi\in H(\operatorname{div};\Omega)$, $j=1,\ldots,d$, where ${\bf e}_{j}$ is the $j$-th unit basis vector of $\mathbb{R}^{d}$. Integration of (3.1b) over $[0,T]$, definition (3.9) with $\ell=0$ and the periodicity of ${\bf q}$ then yield

[TABLE]

From (3.12), we conclude that $\partial_{x_{j}}\gamma_{0}^{v}=0$, $j=1,\ldots,d$, which implies

[TABLE]

Since $\gamma_{0}^{v}$ satisfies (3.1e) with $\ell=0$, $\eta=\gamma_{0}^{v}=0$, if $|\Gamma_{D}|>0$. Similarly, if $|\Gamma_{S}>0|$, (3.1c), together with $\boldsymbol{\gamma}_{0}^{p}\cdot{\bf n}=0$ on $\Gamma_{S}$, yields

[TABLE]

Thus, $\eta=0$ when $|\Gamma_{D}\cup\Gamma_{S}|>0$, which completes the proof.

∎

For sound-soft scattering problems, where $|\Gamma_{D}|>0$ and $|\Gamma_{S}|>0$, $\eta=0$ and all eigenfunctions $\boldsymbol{\gamma}_{\ell}^{p}$, $\gamma_{\ell}^{v}$, $|\ell|>1$ of (3.8) trivially vanish in (3.6) [21]. Therefore, (3.6)-(3.7) in Theorem 2 with $t=0$ imply that

[TABLE]

From the real part of (2.1) we than conclude that

[TABLE]

3.2 Fundamental frequency filtering for first-order formulation

When the CMCG method is applied to the first-order formulation (3.1), any minimizer of $\widehat{J}({\bf p}_{0},v_{0})=0$ generally consists of spurious perturbations $\eta$, $\boldsymbol{\lambda}$ and eigenfunctions $\boldsymbol{\gamma}_{\ell}^{p}$, $\gamma_{\ell}^{v}$ superimposed on the desired (unique) solution $u$ of (2.1). To extract $u$ from $({\bf p}_{0},v_{0})$, we apply a filtering approach, similar to that in Section 2.2, and thereby restore uniqueness. Again, we multiply the Fourier series representation in (3.6) of $v$ by $\operatorname{e}^{i\omega t}$ and integrate over $(0,T)$. Since $\eta$ and $\boldsymbol{\lambda}$ are independent of time, while $\operatorname{e}^{i\omega t}$ is orthogonal to $\operatorname{e}^{i\omega\ell t}$, $|\ell|>1$, all spurious modes vanish and the resulting expression simplifies to:

[TABLE]

which immediately yields

[TABLE]

We summarize this result in the following proposition.

Proposition 2**.**

Let $u\in H^{1}(\Omega)$ be the unique solution of (2.1) and $({\bf p},v)\in\mathcal{Q}$ be a $T$-time periodic solution of (3.1). Then $u$ is given by (3.14) .

3.3 Hybrid DG FE-Discretization

In [12], Glowinski et al. used standard Raviart-Thomas (RT) finite elements to discretize (3.1). Since no mass-lumping is available for RT elements on triangles or tetrahedra [22], each time-step then requires the inversion of the mass-matrix. To avoid that extra computational cost, which strongly impedes parallelization, we instead consider the recent hybrid discontinuous Galerkin (HDG) FEM [13] to discretize (2.1) in its corresponding first-order formulation together with (3.1). Then, the mass-matrix is block-diagonal, with (small and constant) block size equal to the number of dof’s per element, so that the time-stepping scheme becomes truly explicit and inherently parallel.

Let $\mathcal{T}_{h}$ denote a regular triangulation of $\Omega_{h}$, $\mathcal{E}_{h}$ the set of all faces and $\mathcal{P}^{r}$ the space of polynomials of degree $r$. In addition, we define

[TABLE]

[TABLE]

For the time integration of (3.18), we use the standard explicit fourth-order Runge-Kutta (RK4) method.

3.4 Convergence and superconvergence

For a FE discretization with piecewise polynomials of degree $r$, we usually expect convergence as $\mathcal{O}(h^{r+1})$ with respect to the $L^{2}$-norm. For the above HDG discretization, however, an extra power in $h$ can be achieved by applying a cheap local post-processing step [13]. The same (super-) convergence in space of order $r+2$ using only $P^{r}$-FE can be achieved with the CMCG method by applying the post-processing step to the numerical solutions $({\bf p}_{h}^{n_{T}},v_{h}^{n_{T}})$ of (3.1) at the final time $T=n_{T}\Delta t$.

[TABLE]

for any element $K\in\mathcal{T}_{h}$. The new approximate solution $u$ is then given by (3.14) with ${\bf p}$ and $v$ replaced by ${\bf p}_{h}^{n_{T},*}$ and $v_{h}^{n_{T},*}$.

To illustrate the accuracy and verify the expected convergence rates for the various FE discretizations in the CMCG method, we now consider the following one-dimensional solution

[TABLE]

of (2.1) in $\Omega=(0,1)$ with $c=1$, $k=5\pi/4$, $\Gamma_{D}=\{0\}$ and $\Gamma_{S}=\{1\}$. Figure 1 shows the error $\|u-u_{h}\|$ obtained with the CMCG method for the first-order formulation (2.2) and a $P^{2}$-HDG discretization on a sequence of increasingly finer meshes $h=2^{-i}$, $i=3,\ldots,6$. Clearly as we refine the mesh, we always reduce the time-step in the RK4 method to satisfy the CFL stability condition. The CG iteration stops once the tolerance $tol=10^{-12}$ is reached. We also compare the solutions obtained with the CMCG method applied to the second-order formulation using a (continuous) $\mathcal{P}^{2}$ or $\mathcal{P}^{3}$-FEM. All numerical solutions display the expected optimal convergence of order $r+1$ with polynomials of degree $r$, while the first-order HDG approach even achieves superconvergence of order $r+2$, once local post-processing is applied to the final CG iterate.

3.5 Physically bounded domain

In the absence of Dirichlet or impedance boundary conditions, the first-order formulation does not yield the correct minimizer of $J$. As a simple remedy, we proposed in Section 3.2 a filtering procedure which removes the unwanted spurious modes. To illustrate the effectiveness of the filtering procedure, we now consider the exact solution of (2.1)

[TABLE]

in $\Omega=(0,1)$ with homogeneous Neumann boundary conditions and $k=\omega=\pi/4$, $c=1$. Note that $k^{2}$ is not an eigenvalue of (2.6) and therefore the solution of (2.1) is well-posed. However, as $(4k)^{2}=\pi^{2}$ indeed corresponds to the first eigenvalue of the negative Laplacian, the CMCG method in general will not yield the correct (unique) solution – see Theorems 1 and 2. Indeed as shown in Figure 2, the original CMCG method [8] applied to the second-order formulation with the energy functional $J$ in (2.4) does not yield the exact solution of (2.1), unlike the numerical solutions obtained after filtering – see Sections 2.2 and 3.2.

4 Numerical results

Here we present a series of numerical examples that illustrate the accuracy, convergence behavior and parallel performance of the CMCG method. First, we verify that the numerical solution $u_{h}$ of (2.1) obtained with the CMCG method converges to the numerical solution $u_{h}^{*}$ obtained with a direct solver for the same spatial FE discretization as the time step $\Delta t\rightarrow 0$ in the numerical integration of (2.2). Next, we evaluate different stopping criteria for the CG iteration in the CMCG Algorithm from Section 2.3. We also compare the CMCG Algorithm to a long-time solution of the wave equation without controllability (“do-nothing” approach) to demonstrate its effectiveness, in particular for nonconvex obstacles. Moreover, we show how an initial run-up yields a judicious initial guess $(v_{0},v_{1})$ for the CG iteration thereby further accelerating convergence. Finally, we apply the CMCG method to large scale scattering problems on a massively parallel architecture, where the elliptic problem (2.16) is solved in parallel with domain decomposition methods.

4.1 Semi-discrete convergence

First, we consider a simple 1D example to show for a fixed FE-mesh that the numerical solution $u_{h}$, obtained with the CMCG method, converges to the numerical solution $u_{h}^{*}$, obtained with a direct solver, as $\Delta t\rightarrow 0$. Hence we consider the following solution $u$ of (2.1) in $\Omega=(0,1)$ with $\omega=k=6\pi$, $c=1$ and $f\equiv 0$:

[TABLE]

Now, let $u_{h}^{*}(x)$ be the FE Galerkin solution corresponding to the direct solution of the linear system

[TABLE]

resulting from the same standard $H^{1}$-conforming or HDG $P^{2}$-FE discretization of the Helmholtz equation (2.1) in second- or first-order formulation, respectively. For the time integration of (2.2) or (3.1) in the CMCG Algorithm, we use the standard explicit fourth order Runge-Kutta (RK4) method.

Usually we avoid inverting the mass-matrix at each time step via order preserving mass-lumping [23] which, however, introduces an additional spatial discretization error. Here to ensure a consistent comparison, we thus compute $u_{h}$ and $u_{h}^{*}$ both either with, or without, mass-lumping (ML). For the CG iteration, we always choose $v_{0}^{(0)}\equiv 0$, $v_{1}^{(0)}\equiv 0$ and fix the tolerance to $tol=10^{-14}$ to ensure convergence to machine precision accuracy.

In Figure 3, we monitor the difference between the numerical solution $u_{h}^{*}$ or $u_{h,HDG}^{*}$ of (4.1), obtained with a direct solver, and $u_{h}$ or $u_{h,HDG}$, obtained with the CMCG method using either the second or the first order formulation, respectively. As expected, for increasingly smaller $\Delta t$ and a fixed stringent tolerance in the CG iteration, the numerical solution of the CMCG method always converges to the discrete solution of the Helmholtz equation for the same FE discretization.

4.2 CG iteration and initial run-up

Next, we first compare different stopping criteria for the CG iteration in the CMCG Algorithm applied to the original second-order formulation from Section 2. We then illustrate how the CMCG method greatly accelerates the convergence of a solution of the wave equation to its long-time asymptotic limit, in particular for nonconvex obstacles. Finally, we show how an initial run-up yields a judicious initial guess for the CG iteration, which further accelerates the convergence of the CMCG Algorithm.

Hence, we consider a two-dimensional sound-soft scattering problem (2.1) with $c\equiv 1$, $k=\omega=2\pi$, $f\equiv g_{D}\equiv g_{N}\equiv 0$ and $g_{S}=-(\partial_{n}-ik)u^{in}$ in a bounded square domain $\Omega=(0,10\lambda)\times(0,10\lambda)$, $\lambda=1$, either with a convex obstacle or a semi-open square shaped cavity. On the boundary $\Gamma_{D}$ of the obstacle, we impose a homogeneous Dirichlet condition and on the exterior boundary $\Gamma_{S}$ a Sommerfeld-like absorbing condition on the total wave field. The incident plane wave

[TABLE]

impinges with the angle $\theta=135^{\circ}$ upon the obstacle.

4.2.1 CG iteration and stopping criteria

In Algorithm (Section 2.3), the CMCG method terminates at the $\ell$-th iteration and returns

[TABLE]

when the relative CG-residual in Step 5.8,

[TABLE]

is less than the tolerance $tol$. Indeed, a small CG-residual indicates that the gradient of $J$ is sufficiently small at $(v_{0}^{(\ell)},v_{1}^{(\ell)})$ and thus that a minimum has been reached.

Since the cost functional $J$ also vanishes at the minimum, we can use $J$ itself, instead of its gradient, to monitor convergence of the CG iteration via the relative periodicity misfit,

[TABLE]

In fact, the convergence criterion (4.5) is typically used in long-time simulations of the wave equation without controllability (“do-nothing” approach) to determine the current misfit from periodicity in the energy norm.

Alternatively, we may also directly compute the current relative Helmholtz residual from (2.1):

[TABLE]

where ${\bf A}_{h}$ and ${\bf b}_{h}$ result from a FE discretization of (2.1) without mass-lumping, ${\bf u}_{h}^{(\ell)}$ corresponds to the discrete vector of FE coefficients of $u_{h}^{(\ell)}$, and $\|\cdot\|_{2}$ denotes the discrete Euclidean norm.

In Figure 5, we monitor $|u_{h}^{(\ell)}|_{CG}$, $|u_{h}^{(\ell)}|_{J}$ and $|u_{h}^{(\ell)}|_{H}$, defined in (4.4)–(4.6) for the CMCG solution $u_{h}^{(\ell)}$ at the $\ell$-th CG iteration. Whether for a convex (Figure 4a) or a nonconvex (Figure 4b) obstacle, both the CG-residual $|u_{h}^{(\ell)}|_{CG}$ and the periodicity misfit $|u_{h}^{(\ell)}|_{J}$ rapidly converge to zero. In contrast, the Helmholtz residual $|u_{h}^{(\ell)}|_{H}$ stagnates beyond the first hundred CG iterations, as the mass-matrix that appears in ${\bf A}_{h}$ in (4.6) is discretized here without mass-lumping. That additional discretization error together with the numerical error in the time integration of (2.2) both prevent the discrete Helmholtz residual $|u_{h}^{(\ell)}|_{H}$ from converging to zero; hence, (4.6) is generally not a reliable stopping criterion for the CMCG method, unless the spatial FE discretizations used in (2.1) and (4.1) are identical.

4.2.2 CMCG method vs. long-time wave equation solver

In general, the solution $w(x,t)$ of the time-harmonically forced wave equation (2.2) converges asymptotically to the time-harmonic solution [24]

[TABLE]

where $u$ is the (unique) solution of the Helmholtz equation (2.1). Thus, with a wave equation solver at hand, one can in principle compute $u$ from $w$ by solving (2.2) without controllability until a quasi-periodic regime is reached. Given the current value of $w(\cdot,t)$ at time $t=\ell\,T$, $\ell\geq 1$, one can extract from it the complex-valued approximate solution of (2.1),

[TABLE]

which converges to $u$ as $\ell\rightarrow+\infty$. This “do-nothing” approach only requires the time integration of (2.2) without controllability or CG iteration, but it may converge arbitrarily slowly for nonconvex obstacles due to trapped modes [8, 11].

In Figure 6, we monitor the periodicity misfit of $|u_{h}^{(\ell)}|_{J}$ and $|w_{h}^{(\ell)}|_{J}$, where $u_{h}^{(\ell)}$ is the CMCG solution at the $\ell$-th CG iteration and $w_{h}^{(\ell)}$ is given by (4.8). In addition, we also compare both numerical solutions with the direct solution $u_{h}^{*}$ of the linear system (4.1), resulting from the same underlying FE discretization, yet without mass-lumping.

We observe that the asymptotic solution $w_{h}^{(\ell)}$ and the CMCG solution $u_{h}^{(\ell)}$ indeed both converge to the time-harmonic solution $u_{h}^{*}$, until the additional errors caused by mass-lumping and the time discretization dominate the total error – see Section 4.1. For the convex obstacle, the number of CG iterations required by $u_{h}^{(\ell)}$ is only half the number of time periods needed for $w_{h}^{(\ell)}$ to reach the same level of accuracy. However, since each CG iteration requires not only the solution of a forward and backward wave equation but also of the elliptic problem (2.16a), simply computing a long-time solution of the time-harmonically forced wave equation (2.2) without controllability in fact proves cheaper here than the CMCG Algorithm. For a nonconvex obstacle, however, the long-time numerical solution of the time-dependent wave equation $w_{h}^{(\ell)}$ converges extremely slowly and fails to reach the asymptotic time-harmonic regime even after 1000 periods. In contrast, the convergence of the CMCG solution $u_{h}^{(\ell)}$ remains remarkably insensitive to the non-convexity of the obstacle.

4.2.3 Initial run-up

In [25], Mur suggested that convergence of the time-harmonically forced wave equation (2.2) to the time-harmonic asymptotic regime can be accelerated by pre-multiplying the time-harmonic sources in (2.2) with the smooth transient function $\theta_{tr}$ from zero to one,

[TABLE]

active during the initial time interval $[0,t_{tr}]$, $t_{tr}=\ell\,T$ – see also [8].

Again, we consider plane wave scattering either from a convex or nonconvex obstacle – see Figure 4. Now, we first solve the wave equation (2.2) with the modified source terms and zero initial conditions until time $t=\ell\,T$, $\ell\geq 1$, which yields the time-dependent solution $y_{tr}$. After that initial run-up phase, we then apply the CMCG Algorithm (Section 2.3) using the initial guess

[TABLE]

To estimate the total computational effort, we count the total number of time periods for which the (forward or backward) wave equation is solved: $\ell$ during initial run-up and $2\times\#iter_{CG}$ during the CG iteration. In Figure 7 we display the total number $2\times\#iter_{CG}+\ell$ of time periods needed until convergence with $tol=10^{-6}$, as we vary the number of periods $\ell$ in the initial run-up.

For a convex obstacle, the CMCG Algorithm without any initial run-up requires $888$ time periods. However, as in Section 4.2, convergence can also be achieved at a comparable computational effort simply by solving the wave equation, here with the source terms pre-multiplied by $\theta_{tr}$ in (4.9). Still, the minimal computational cost is achieved when both the initial run-up and the CMCG Algorithm are combined.

For the nonconvex obstacle, however, simply solving the time-harmonically forced wave equation over a very long time, be it with or without $\theta_{tr}(t)$ smoothing, fails to reach the long-time asymptotic final time-harmonic state. Regardless of the length of the initial run-up, convergence indeed cannot be achieved here (within 1000 time periods) without controllability because of trapped modes. Nevertheless, the initial run-up always speeds up the convergence of the CMCG method by providing a judicious initial guess for the CG iteration.

4.3 Parallel computations

Both the CMCG method for the second-order formulation from Section 2 and that for the first-order formulation from Section 3 lead to inherently parallel non-intrusive algorithms, as long as an efficient parallel solver for the time-dependent wave equation is available. As the first-order formulation with the HDG discretization neither requires mass-lumping nor the solution of an elliptic problem, it is in fact trivially parallel. Here we demonstrate that even the CMCG approach for the second-order formulation, which does require the solution of (2.16a) at each CG iteration, nonetheless achieves strong scalability on a massively parallel architecture.

The CMCG Algorithm from Section 2.3 is implemented within FreeFem++ [26], an open source finite element software written in C++. FreeFem++ defines a high-level Domain Specific Language (DSL) and natively supports distributed parallelism with MPI. The parallel implementation of the CMCG method relies on the spatial decomposition of the computational domain $\Omega$ into multiple subdomains, each assigned to a single computing core. Local finite element spaces are then defined on the local meshes of the subdomains, effectively distributing the global set of degrees of freedom across the available cores.

The bulk of the computational work for solving the forward and backward wave equations in Step (5)5.1  of the CMCG Algorithm simply consists in performing a sparse matrix-vector product at each time step, which is easily parallelized in this domain decomposition framework: it amounts to performing local matrix-vector products in parallel on the local set of degrees of freedom corresponding to each subdomain, followed by local exchange of shared values between neighboring subdomains.

While the explicit time integration of the wave equation is trivially parallelized thanks to mass-lumping, achieving good parallel scalability for the elliptic problem in Step (5)5.2  of the CMCG Algorithm is more difficult. Here we use domain decomposition (DD) methods [18], which are well-known to produce robust and scalable parallel preconditioners for the iterative solution of large scale partial differential equations. We use the parallel DD library HPDDM [27], which implements efficiently various Schwarz and substructuring methods in C++11 with MPI and OpenMP for parallelism and is interfaced with FreeFem++ .

The elliptic problem (2.16a) in the CMCG algorithm is solved by HPDDM using a two-level overlapping Schwarz DD preconditioner, where the coarse space is built using Generalized Eigenproblems in the Overlap (GenEO) [28]. The GenEO approach has proved effective in producing highly scalable preconditioners for solving various elliptic problems  [6, 28].

All computations were performed on the supercomputer OCCIGEN at CINES, France 111https://www.cines.fr/calcul/materiels/occigen/, with 50544 (Intel XEON Haswell) cores.

4.3.1 2D Marmousi Model

Here we consider the well-known Marmousi model from geophysics [29], that is (2.1) in $\Omega=(0,9.2)\times(0,3)$ $[km]$ with the source

[TABLE]

The velocity profile $c(x)$ is shown in Figure 8 and we apply absorbing boundary conditions on the lateral and lower boundaries and a homogeneous Dirichlet condition at the top. For the spatial discretization, we use a $P^{2}$-FE method with (order preserving) mass-lumping [23] and at least $15$ points per wave length. For the time integration of (2.2), we apply the leap-frog scheme (LF); here, the number of $T/\Delta t=390$ time steps per period remains constant at all frequencies $\nu=\omega/2\pi$, as both $T$ and $\Delta t$ are inversely proportional to $\nu$. To speed-up the convergence of the CMCG method, we also use an initial run-up (Section 4.2) until time $t_{tr}$, which lets waves travel at least once across the entire computational domain during run-up; hence, we set

[TABLE]

For any particular frequency $\nu$, we apply the CMCG method for fixed parameters and FE-mesh while increasing the number of (CPU) cores. Figure 9 displays the real part of the wave field with $\nu=250$ [Hz]. In Figure 10, we observe linear speed-up (strong scaling) at every frequency with increasing number of cores. In fact, the speed-up is even slightly better than linear due to cache effects, but also because the cost of the direct solver used on each subdomain decreases superlinearly with the decreasing size of subdomains as the number of cores increases.

As the frequency $\nu$ increases, both the period $T=1/\nu$ and the time-step $\Delta t$ decrease, so that the number of time steps per CG iteration remains constant. Since the number of CG iterations does not grow here with increasing $\nu$, the bulk of the computational work in the CMCG Algorithm in fact shifts to the run-up phase. For $\nu=10$ Hz, for instance, the CMCG Algorithm stops after $273$ CG iterations, while $74\%$ of the total computational time is spent in the time integration of (2.2), $16\%$ in the elliptic solver (DDM) and $10\%$ in the initial run-up. In contrast, for $\nu=250$ Hz, the CMCG Algorithm already stops after $5$ CG iterations, while $99\%$ of the total computational time is spent in the initial run-up and $1\%$ in the CG iteration. By modifying the run-up time $t_{tr}$, one could arbitrarily shift the relative computational cost between run-up and CG iterations and thus further optimize for a minimal total execution time.

4.3.2 3D cavity

Finally, we compute the scattered wave from a sound-soft cavity – see Figure 11 – and hence consider (2.1) in $\Omega=(0,6)\times(0,3)\times(0,3)$ with $c=1$, $k=\omega=2\pi\nu$, $\lambda=1$, $f\equiv g_{D}\equiv g_{N}\equiv 0$ and

[TABLE]

We impose a homogeneous Dirichlet boundary condition on the obstacle and a Sommerfeld-like absorbing condition on the exterior boundary of $\Omega$.

Now, we discretize (2.2) with $P^{1}$-FE in space and the second-order LF method in time. To control the pollution error, we set $hk^{3/2}\sim const$, as we increase the frequency $\nu$. Figure 12 shows the total wave field with $\nu=6$ inside the cavity. For fixed parameters and mesh size, we now solve (2.1) at frequencies $\nu=2,3,4,6$ with the CMCG method using an increasing number of cores – see Table 2. Again, we observe in Figure 13 (better than) linear (strong) scaling with increasing number of cores. In contrast to the previous Marmousi problem, the ”do-nothing” approach without controllability fails here because the 3D cavity is not convex.

5 Concluding remarks

We have presented two inherently parallel controllability methods (CM) for the numerical solution of the Helmholtz equation in heterogeneous media. The first, based on the second-order formulation of the wave equation, uses a standard (continuous) FE discretization in space with order preserving mass-lumping. Each conjugate gradient (CG) iteration then requires the explicit time integration of a forward and backward wave equation, together with the solution of the symmetric and coercive elliptic problem (2.16), which is independent of the frequency. The second, based on the first-order (or mixed) formulation of the wave equation, uses a recent hybridized discontinuous Galerkin (HDG) discretization, which not only automatically yields a block-diagonal mass-matrix but also completely avoids solving (2.16). Hence, it is trivially parallelized and even leads to superconvergence after a local post-processing step.

Both CMCG methods are inherently parallel, as they lead to iterative algorithms whose convergence rate is independent of the number of cores on a distributed memory architecture. Thanks to the well-known parallel efficiency of explicit methods combined with the excellent scalability of two-level domain decomposition preconditioners for coercive elliptic problems up to thousands of cores implemented in HPDDM, even the second-order CMCG approach exhibits parallel strong scalability.

The CMCG method can be applied to general boundary-value problems governed by the Helmholtz equation, such as sound-soft or sound-hard scattering problems or wave propagation in physically bounded domains. Although the CMCG solution will generally contain higher order spurious eigenmodes, we have proposed in Section 2.2 a simple filtering procedure to remove them. Furthermore, including a transient initial run-up to determine a judicious initial guess significantly accelerates the CG iteration. In fact, for scattering from convex obstacles, simply solving the time-harmonically forced wave equation over a long-time without any controllability can provide an even simpler, highly parallel Helmholtz solver. For nonconvex obstacles, however, solving the wave equation without any controllability (”do-nothing” approach) is not a viable option, as the long time asymptotic convergence to the time-harmonic regime is simply too slow due to trapped modes. In all cases, the CMCG Algorithm combined with the initial run-up leads to the smallest time-to-solution.

The CMCG approach developed here for the Helmholtz equation immediately generalizes to other time-harmonic vector wave equations from electromagnetics or elasticity. Its implementation is non-intrusive and particularly useful when a parallel efficient time-dependent wave equation solver is at hand. In the presence of local mesh refinement, local time-stepping methods [31] permit to circumvent the increasingly stringent CFL condition without sacrificing the explicitness or inherent parallelism. Finally, the CMCG method can also be used to compute periodic, but not necessarily time-harmonic, solutions of the wave equations. In particular, if the source consists of a superposition of several time-harmonic sources (”super-shot”) with rational frequencies, the solutions to the different Helmholtz problems can be extracted via filtering from a single application of the CMCG method.

Acknowledgement: This work was supported by the Swiss National Science Foundation under grant SNF 200021_169243. Access to the HPC resources of CINES was granted under allocation 2018-A0040607330 by GENCI.