Parallel Controllability Methods
For the Helmholtz Equation
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
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
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
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
To determine and , the problem is reformulated as a least-squares optimization problem over for the quadratic cost functional
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 .
for any , where the constants and the eigenfunctions , , satisfy
Let . Then satisfies
Furthermore, if , then . If , then .
Here and denotes the standard inner product.
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
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
where and all have vanished but the constant is still undetermined.
from (2.1a). This immediately yields the remaining constant
We summarize the above derivation in the following proposition.
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 :
Here denotes an arbitrary perturbation, the standard duality pairing, and the solution of the adjoint (backward) wave equation:
and the initial conditions satisfy for any
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:
For the sake of completeness, we list the full CMCG Algorithm – see BGP1998; GT2018:
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:
with the initial conditions
Hence, the solution of (3.1) lies in the function space E2010; P1983,
In terms of and , the energy functional defined in (2.2) now becomes
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 :
Here denotes an arbitrary perturbation with
whereas solves the backward (adjoint) wave equation in first-order form GR2006, that is (3.1) with and
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.
where the constant , with
and the complex-valued eigenfunctions , , satisfy
Furthermore, if , then .