Parallel Controllability Methods For the Helmholtz Equation

# Parallel Controllability Methods For the Helmholtz Equation

Marcus J. Grote Frédéric Nataf Jet Hoe Tang Pierre-Henri Tournier University of Basel, Spiegelgasse 1, 4051 Basel, Switzerland Laboratoire J.L. Lions, Université Pierre et Marie Curie, 4 place Jussieu, 75005 Paris, France, and ALPINES INRIA, Paris, France
###### 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.

###### keywords:
Helmholtz equation; time-harmonic scattering; exact controllability; finite elements; domain decomposition; parallel scalability
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 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 (EG2005).

In recent years, a growing number of increasingly sophisticated preconditioners has been proposed for the iterative solution of the Helmholtz equation; ”Shifted Laplacian” preconditioners ERLANGGA2004409, for instance, have led to modern multigrid Calandra:2017:GMP; BGS2009 and domain decomposition preconditioners graham2017domain; BDGST2017. While some of those preconditioners may achieve a desirable frequency independent convergence behavior in special situations EY2011, 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, , with known period . By minimizing an energy functional which penalizes the mismatch after one period, controllability methods iteratively adjust the (unknown) initial condition thereby steering towards the desired periodic solution. Once the minimizer of 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 BGP1998, 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 then lies in , the original formulation requires the solution of a coercive elliptic problem at each CG iteration. Heikkola et al. in HMPR2007; HMPR2007_2 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 is no longer unique. In GT2018, 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 GR2006 proposed a CM based on the wave equation in first-order (or mixed) form using classical Raviart-Thomas (RT) finite elements. As then lies in , 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 BGP1998; GR2006 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 , 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 (CNPS2016), 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 in a bounded connected computational domain , , with a Lipschitz boundary . The boundary consists of three disjoint components, 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 , the wave field hence satisfies the Helmholtz equation

 −Δu(x) −k2(x) u(x) = f(x), x∈Ω, (2.1a) ∂u(x)∂n − ik(x) u(x) = gS(x), x∈ΓS, (2.1b) ∂u(x)∂n = gN(x), x∈ΓN, (2.1c) u(x) = gD(x), x∈ΓD, (2.1d)

where is the (angular) frequency, the wave speed, the wave number, the unit outward normal, and , , and are known and may vanish.

The above formulation is rather general and encompasses most common applications such as sound-soft scattering problems with and , sound-hard scattering problems with and , or physically bounded domains with . We shall always assume for any particular choice of , , or combination of boundary conditions that (2.1) is well-posed and has a unique solution .

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, , satisfies the (real-valued) time-dependent wave equation

 1c2(x)∂2y(x,t)∂2t − Δy(x,t) = Re{f(x)\e−iωt}, x∈Ω, t>0, (2.1a) ∂y(x,t)∂n+1c(x) ∂y(x,t)∂t = Re{gS(x)\e−iωt}, x∈ΓS,t>0, (2.1b) ∂y(x,t)∂n = Re{gN(x)\e−iωt}, x∈ΓN,t>0, (2.1c) y(x,t) = Re{gD(x)\e−iωt}, x∈ΓD,t>0, (2.1d) y(x,0) = v0(x),∂y(x,0)∂t = v1(x), x∈Ω, (2.1e)

for the (unknown) initial values and .

For sound-soft scattering problems (2.1), where and , Bristeau et al. BGP1998; L1988 proposed to determine via controllability by computing a time-periodic solution of (2.1) with period . Once the initial values of are known, the solution of the original Helmholtz equation (2.1) is immediately given by

 u=v0+iωv1,v0,v1∈H1(Ω). (2.1)

To determine and , the problem is reformulated as a least-squares optimization problem over for the quadratic cost functional

 J(v0,v1)=12∫Ω|∇y(x,T)−∇v0(x)|2dx+12∫Ω1c2(x)(yt(x,T)−v1(x))2dx, (2.2)

where satisfies (2.1) with the initial values and . The functional measures in the energy norm the mismatch between the solution of (2.1) at the initial time and after one period. It is non-negative and convex, while if, and only if, and for any given initial values ; in particular, if and .

For more general scattering problems, however, is no longer strictly convex as the -periodicity of and no longer guarantees a unique periodic solution of (2.1). Instead, for the general boundary-value problem (2.1), the situation is more complicated and summarized in the following theorem GT2018 .

###### Theorem 1.

Let be the unique solution of (2.1) and be a (real-valued) solution of (2.1) with initial values . If and are time periodic with period , then admits the Fourier series representation

 (y(⋅,t),φ) = (Re{u\e−iωt},φ)+(λ+ηt,φ)+∑|ℓ|>1(γℓ\eiωℓt,φ) (2.3)

for any , where the constants and the eigenfunctions , , satisfy

 −Δγℓ(x)−(ℓk(x))2γℓ(x) = 0, x∈Ω, (2.4a) ∂γℓ(x)∂n+iℓk(x) γℓ(x) = 0, x∈ΓS, (2.4b) ∂γℓ(x)∂n = 0, x∈ΓN, (2.4c) γℓ(x) = 0, x∈ΓD, (2.4d)

Let . Then satisfies

 (v,φ) = (u,φ)+(λ+iωη,φ)+∑|ℓ|>1(αℓ+iℓβℓ,φ),∀φ∈H1D. (2.1)

Furthermore, if , then . If , then .
Here and denotes the standard inner product.

###### Proof.

See GT2018. ∎

For sound-soft scattering problems (), where both Dirichlet and Sommerfeld-like absorbing boundary conditions are imposed on , all the eigenfunctions , , and the constants , in (2.1) vanish identically. Thus, the minimizer of in (2.2) then coincides with .

For scattering problems from sound-hard obstacles or penetrable inclusions (, ), the eigenfunctions and the constant in (2.1) still vanish identically, yet the constant may be nonzero. Given any minimizer of , we can recover by subtracting the spurious shift using the compatibility condition:

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 BGT1982; GK1995, also circumvents the indeterminacy due to . For wave propagation in physically bounded domains (), the eigenfunctions and the constants in (2.1) typically do not vanish. However, we can always restore uniqueness by replacing with an alternative energy functional, thereby incurring a small increase in computational cost – see GT2018.

### 2.2 Fundamental frequency extraction via filtering

From Theorem 1 we conclude that a minimizer of generally yields a time-dependent solution of (2.1), which contains a constant shift determined by , a linearly growing part determined by , and higher frequency harmonics determined by , all superimposed on the desired time-harmonic field with fundamental frequency . Those spurious modes can be eliminated by replacing with an alternative energy functional at a small extra computational cost GT2018. Instead we now propose an alternative approach via filtering which removes all spurious modes without requiring a modified energy functional.

Let be the time-dependent solution of (2.1) that corresponds to a minimizer of . Next, we define as

 ˆy(x):=1T∫T0(y(x,t)+iωyt(x,t))\eiωtdt. (2.2)

To extract from , we now take advantage of the mutual orthogonality of different time harmonics in . Hence, we multiply (2.3) with and integrate in time over to obtain

 ˆy(x) = 1T∫T0(Re{u\e−iωt}+λ+ηt+iIm{u\e−iωt}+iηω)\eiωtdt = 1T∫T0u\e−iωt\eiωtdt−iηω = u−iηω. (2.3)

This yields

 u(x)=ˆy(x)+iηω,x∈Ω (2.4)

where and all have vanished but the constant is still undetermined.

If or , Theorem 1 implies that and thus . Otherwise in the pure Neumann case (), we determine by integrating (2.4), multiplied by , over and using the compatibility condition

 −∫Ωk2(x)u(x) dx=∫Ωf(x) dx+∫∂ΩgN(x) ds. (2.5)

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

 iηω = −1∥k∥2L2(Ω)(∫Ωf(x) dx+∫∂ΩgN(x) ds+∫Ωk2(x)ˆy(x)dx). (2.6)

We summarize the above derivation in the following proposition.

###### Proposition 1.

Let be the unique solution of (2.1) and the time dependent solution of (2.1) corresponding to a minimizer of , i.e. . Then is given by (2.4) with if or , and with given by (2.6) when .

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

### 2.3 The CMCG Algorithm

To minimize the quadratic cost functional defined by (2.2) over , a natural choice is the conjugate gradient (CG) method BGP1998, which requires the Fréchet derivative of at :

 ⟨J′(v),δv⟩ = −∫Ω∇(y(x,T)−v0(x))⋅∇δv0(x) dx −∫Ω1c2(x)(yt(x,T)−v1(x))δv1(x) dx +∫Ω1c2(x)(p(x,0)δv1(x)−pt(x,0)δv0(x)) dx +∫ΓS1c(x) p(x,0)δv0(x) ds.

Here denotes an arbitrary perturbation, the standard duality pairing, and the solution of the adjoint (backward) wave equation:

 1c2(x)∂2∂2tp(x,t) − Δp(x,t) = 0, x∈Ω, t>0, (2.9a) ∂p(x,t)∂n−1c ∂∂tp(x,t) = 0, x∈ΓS,t>0 (2.9b) ∂p(x,t)∂n = 0, x∈ΓN,t>0, (2.9c) p(x,t) = 0, x∈ΓD,t>0, (2.9d) p(x,T) = p0(x),∂p(x,T)∂t = p1(x), x∈Ω, (2.9e)

and the initial conditions satisfy for any

 p0(x) = yt(x,T)−v1(x),x∈Ω, ∫Ωp1(x)c2(x)w(x) dx = ∫ΓSp0(x)c(x)w(x) ds−∫Ω∇(y(x,T)−v0(x))⋅∇w(x) dx.

The derivation of (2.3) and (2.3) can be found in BGP1998. In each CG iteration the derivative requires the solution of the forward and backward (adjoint) wave equations (2.1) and (2.3) over one period . Moreover, each CG iteration requires an explicit (Riesz) representer of the gradient defined in (2.3), which is determined by solving the symmetric and coercive elliptic problem BGP1998; MS2014:

 (∇~g0,∇φ) = ∫Ωg0(x)φ(x) dx = ∫Ω∇(v0(x)−y(x,T))⋅∇φ(x)−1c2(x)pt(x,0)φ(x) dx +∫ΓS1c(x) p(x,0)φ(x) ds,∀φ∈H1D, (2.1a) ~g1 = g1 = v1−yt(⋅,T)+c−2p(⋅,0). (2.1b)

For the sake of completeness, we list the full CMCG Algorithm – see BGP1998; GT2018:

CMCG Algorithm.

1. Initialize (initial guess).

2. Solve the forward and the backward wave equations (2.1) and (2.3) to determine the gradient of , , defined by (2.3).

3. Solve the coercive elliptic problem (2.3) with to determine the new search direction .

4. Set .

5. For

1. Solve the wave equation (2.1) with and the initial values and the backward wave equation (2.3). Compute the gradient defined by (2.3).

2. Solve the coercive elliptic problem (2.3) with to get .

3. Stop when the relative residual lies below the given tolerance

    ⎷∥∇r(ℓ+1)0∥2L2(Ω)+∥(1/c) r(ℓ+1)1∥2L2(Ω)∥∇r(0)0∥2L2(Ω)+∥(1/c) r(0)1∥2L2(Ω)≤tol.
6. Return approximate solution of (2.1) given by

 uh=v(ℓ)0+iωv(ℓ)1.

Since , the updates of , and in Steps 5.4, 5.5 and 5.7 in the CMCG Algorithm also remain in . We emphasize that (2.1a) is independent of 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 DJN2015; BDGST2017.

## 3 Controllability methods for first-order formulations

The CMCG Algorithm from Section 2.3 iterates on the initial value of the second-order wave equation (2.1) until its solution is -time periodic. However, the gradient of the cost functional , which is needed during the CG update, only lies in the dual space . To ensure that the solution remains sufficiently regular and in , the corresponding Riesz representative is computed at every CG iteration by solving the strongly elliptic problem (2.1a). In GR2006, Glowinski et al. derived an equivalent first-order formulation for sound-soft scattering problems, where the solution instead lies in , which is reflexive. As a consequence, all CG iterates automatically lie in the correct solution space , while the solution of (2.1a) is no longer needed.

### 3.1 First-order formulation for general boundary conditions

Again, we always assume for any particular choice of , , and combination of boundary conditions that (2.1) has a unique solution . Following GR2006, we now let , and rewrite the time-dependent wave equation (2.1) in first-order form:

 1c2(x)vt(x,t)−∇⋅p(x,t) = Re{f(x)\e−iωt}, x∈Ω, t>0, (3.1a) ∂∂tp(x,t) = ∇v(x,t), x∈Ω, t>0, (3.1b) p(x,t)⋅n+1c(x)v(x,t) = Re{gS(x)\e−iωt}, x∈ΓS,t>0, (3.1c) p(x,t)⋅n = Re{gN(x)\e−iωt}, x∈ΓN,t>0, (3.1d) v(x,t) = Re{−iωgD(x)\e−iωt}, x∈ΓD,t>0 (3.1e)

with the initial conditions

 p(x,0) = p0(x)∈Rd,v(x,0) = v0(x)∈R, x∈Ω. (3.1f)

Hence, the solution of (3.1) lies in the function space  E2010; P1983,

 Q = C0([0,T];H(div;Ω)×L2(Ω))∩C1([0,T];(L2(Ω))d+1). (3.1)

In terms of and , the energy functional defined in (2.2) now becomes

 ˆJ(p0,v0)=12∫Ω|p(x,T)−p0(x)|2 dx+12∫Ω1c2(x)(v(x,T)−v0(x))2 dx, (3.2)

where solves (3.1) with initial value .

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 is now replaced by :

 ⟨ˆJ′(p0,v0),(δp0,δv0)⟩= ∫Ω(p∗(x,0)−p∗(x,T))δp0(x) dx (3.3) +∫Ω(v∗(x,0)−v∗(x,T))δv0(x) dx.

Here denotes an arbitrary perturbation with

 P={p∈H(div;Ω) | p⋅n=0 on ΓN}, (3.4)

whereas solves the backward (adjoint) wave equation in first-order form GR2006, that is (3.1) with and

 p∗(⋅,T)=p(⋅,T)−p0,v∗(⋅,T)=v(⋅,T)−v0.

For sound-soft scattering problems (), the functional always has a unique (global) minimizer, which therefore coincides with the (unique) time-harmonic solution of (3.1). For more general boundary value problems, however, the minimizer of is not necessarily unique, as shown in the following theorem.

###### Theorem 2.

Let be the unique solution of (2.1) and be a real-valued solution of (3.1) with initial values . If and are time periodic with period , then and admit the Fourier series representation

 p(⋅,t) =Re{∇u\e−iωt}+λ+∞∑|ℓ|>1γpℓ\e−iωℓt, (3.5a) v(⋅,t) =ωIm{u\e−iωt}+η+∞∑|ℓ|>1γvℓ\e−iωℓt, (3.5b)

where the constant , with

 ∫Ωλ⋅∇φ dx=0,∀φ∈H1(Ω),φ|ΓD≡0, (3.1)

and the complex-valued eigenfunctions , , satisfy

 −c2(x)∇⋅γpℓ(x)+iωℓγvℓ(x) = 0, x∈Ω, (3.2a) iωℓγpℓ(x) = ∇γvℓ(x), x∈Ω, (3.2b) c(x)γpℓ(x)⋅n+γvℓ(x) = 0, x∈ΓS, (3.2c) γpℓ(x)⋅n = 0, x∈ΓN, (3.2d) γvℓ(x) = 0, x∈ΓD. (3.2e)

Furthermore, if , then .

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters