11institutetext: 1Department of Applied Mathematics, Naval Postgraduate
School, Monterey, CA
2Institute for Computational Engineering & Sciences, The University of Texas
at Austin, Austin, TX
3Department of Aerospace Engineering & Engineering Mechanics, The University
of Texas at Austin, Austin, TX
4Departments of Mechanical Engineering and Jackson School of Geosciences, The
University of Texas at Austin, Austin, TX

# Discretely exact derivatives for hyperbolic PDE-constrained optimization
problems discretized by the discontinuous Galerkin method111This document has
been approved for public release; its distribution is unlimited.

Lucas C. Wilcox1 Georg Stadler2
Tan Bui-Thanh2,3 and Omar Ghattas2,4

###### Abstract

This paper discusses the computation of derivatives for optimization problems
governed by linear hyperbolic systems of partial differential equations (PDEs)
that are discretized by the discontinuous Galerkin (dG) method. An efficient
and accurate computation of these derivatives is important, for instance, in
inverse problems and optimal control problems. This computation is usually
based on an adjoint PDE system, and the question addressed in this paper is
how the discretization of this adjoint system should relate to the dG
discretization of the hyperbolic state equation. Adjoint-based derivatives can
either be computed before or after discretization; these two options are often
referred to as the optimize-then-discretize and discretize-then-optimize
approaches. We discuss the relation between these two options for dG
discretizations in space and Runge–Kutta time integration. The influence of
different dG formulations and of numerical quadrature is discussed. Discretely
exact discretizations for several hyperbolic optimization problems are
derived, including the advection equation, Maxwell’s equations and the coupled
elastic-acoustic wave equation. We find that the discrete adjoint equation
inherits a natural dG discretization from the discretization of the state
equation and that the expressions for the discretely exact gradient often have
to take into account contributions from element faces. For the coupled
elastic-acoustic wave equation, the correctness and accuracy of our derivative
expressions are illustrated by comparisons with finite difference gradients.
The results show that a straightforward discretization of the continuous
gradient differs from the discretely exact gradient, and thus is not
consistent with the discretized objective. This inconsistency may cause
difficulties in the convergence of gradient based algorithms for solving
optimization problems.

###### Keywords:

Discontinuous Galerkin PDE-constrained optimization Discrete adjoints Elastic
wave equation Maxwell’s equations

##  1 Introduction

Derivatives of functionals, whose evaluation depends on the solution of a
partial differential equation (PDE), are required, for instance, in inverse
problems and optimal control problems, and play a role in error analysis and a
posteriori error estimation. An efficient method to compute derivatives of
functionals that require the solution of a state PDE is through the solution
of an adjoint equation. In general, this adjoint PDE differs from the state
PDE. For instance, for a state equation that involves a first-order time
derivative, the adjoint equation must be solved backwards in time.

If the partial differential equation is hyperbolic, the discontinuous Galerkin
(dG) method is often a good choice to approximate the solution due to its
stability properties, flexibility, accuracy and ease of parallelization.
Having chosen a dG discretization for the state PDE, the question arises how
to discretize the adjoint equation and the expression for the gradient, and
whether and how their discretization should be related to the discretization
of the state equation. One approach is to discretize the adjoint equation and
the gradient independently from the state equation, possibly leading to
inaccurate derivatives as discussed below. A different approach is to derive
the discrete adjoint equation based on the discretized state PDE and a
discretization of the cost functional. Sometimes this latter approach is
called the discretize-then-optimize approach, while the former is known as
optimize-then-discretize. For standard Galerkin discretizations, these two
possibilities usually coincide; however, they can differ, for instance, for
stabilized finite element methods and for shape derivatives
CollisHeinkenschloss02 ; Gunzburger03 ; HinzePinnauUlbrichEtAl09 ; Braack09 .
In this paper, we study the interplay between these issues for derivative
computation in optimization problems and discretization by the discontinuous
Galerkin method.

While computing derivatives through adjoints on the infinite-dimensional level
and then discretizing the resulting expressions (optimize-then-discretize)
seems convenient, this approach can lead to inaccurate gradients that are not
proper derivatives of any optimization problem. This can lead to convergence
problems in optimization algorithms due to inconsistencies between the cost
functional and gradients Gunzburger03 . This inaccuracy is amplified when
inconsistent gradients are used to approximate second derivatives based on
first derivatives, as in quasi-Newton methods such as the BFGS method.
Discretizing the PDE and the cost functional first (discretize-then-optimize),
and then computing the (discrete) derivatives guarantees consistency. However,
when an advanced discretization method is used, computing the discrete
derivatives can be challenging. Thus, understanding the relation between
discretization and adjoint-based derivative computation is important. In this
paper, we compute derivatives based on the discretized equation and then study
how the resulting adjoint discretization relates to the dG discretization of
the state equation, and study the corresponding consistency issues for the
gradient.

_Related work:_ Discretely exact gradients can also be generated via
algorithmic differentiation (AD) GriewankWalther08 . While AD guarantees the
computation of exact discrete gradients, it is usually slower than hand-coded
derivatives. Moreover, applying AD to parallel implementations can be
challenging UtkeHascoetHeimbachEtAl09 . The notion of _adjoint consistency_
for dG (see Hartmann07 ; AlexeSandu10 ; HarrimanGavaghanSuli04 ;
OliverDarmofal09 ; SchutzMay13 ) is related to the discussion in this paper.
Adjoint consistency refers to the fact that the exact solution of the dual (or
adjoint) problem satisfies the discrete adjoint equation. This property is
important for dG discretizations to obtain optimal-order $L^{2}$-convergence
with respect to target functionals. The focus of this paper goes beyond
adjoint consistency to consider consistency of the gradient expressions, and
considers in what sense the discrete gradient is a discretization of the
continuous gradient. For discontinuous Galerkin discretization, the latter
aspect is called dual consistency in AlexeSandu10 . A systematic study
presented in Leykekhman12  compares different dG methods for linear-quadratic
optimal control problems subject to advection-diffusion-reaction equations. In
particular, the author targets commutative dG schemes, i.e., schemes for which
dG discretization and the gradient derivation commute. Error estimates and
numerical experiments illustrate that commutative schemes have desirable
properties for optimal control problems.

_Contributions:_ Using example problems, we illustrate that the discrete
adjoint of a dG discretization is, again, a dG discretization of the
continuous adjoint equation. In particular, an upwind numerical flux for the
hyperbolic state equation turns into a downwind flux in the adjoint, which has
to be solved backwards in time and converges at the same convergence order as
the state equation. We discuss the implications of numerical quadrature and of
the choice of the weak or strong form of the dG discretization on the adjoint
system. In our examples, we illustrate the computation of derivatives with
respect to parameter fields entering in the hyperbolic system either as a
coefficient or as forcing terms. Moreover, we show that discretely exact
gradients often involve contributions at element faces, which are likely to be
neglected in an optimize-then-discretize approach. These contributions are a
consequence of the discontinuous basis functions employed in the dG method and
since they are at the order of the discretization error, they are particularly
important for not fully resolved problems.

_Limitations:_ We restrict ourselves to problems governed by _linear_
hyperbolic systems. This allows for an explicit computation of the upwind
numerical flux in the dG method through the solution of a Riemann problem.
Linear problems usually do not require flux limiting and do not develop shocks
in the solution, which makes the computation of derivatives problematic since
numerical fluxes with limiters are often non-differentiable and defining
adjoints when the state solution involves shocks is a challenge GilesUlbrich10
; GilesUlbrich10a .

_Organization:_ Next, in Section 2, we discuss the interplay of the derivative
computation of a cost functional with the spatial and temporal discretization
of the governing hyperbolic system; moreover, we discuss the effects of
numerical quadrature. For examples of linear hyperbolic systems with
increasing complexity we derive the discrete adjoint systems and gradients in
Section 3, and we summarize important observations. In Section 4, we
numerically verify our expressions for the discretely exact gradient for a
cost functional involving the coupled acoustic-elastic wave equation by
comparing to finite differences, and finally, in Section 5, we summarize our
observations and draw conclusions.

##  2 Cost functionals subject to linear hyperbolic systems

###  2.1 Problem formulation

Let $\Omega\subset\mathbb{R}^{d}$ ($d=1,2,3$) be an open and bounded domain
with boundary $\Gamma=\partial\Omega$, and let $T>0$. We consider the linear
$n$-dimensional hyperbolic system

$\displaystyle\boldsymbol{q}_{t}+\nabla\cdot(\mathbf{F}\boldsymbol{q})$
$\displaystyle=\boldsymbol{f}\qquad$ $\displaystyle\text{\ on\
}\Omega\times(0,T),$ (1a) where, for $(\boldsymbol{x},t)\in\Omega\times(0,T)$,
$\boldsymbol{q}(\boldsymbol{x},t)\in\mathbb{R}^{n}$ is the vector of state
variables and $\boldsymbol{f}(\boldsymbol{x},t)\in\mathbb{R}^{n}$ is an
external force. The flux $\mathbf{F}\boldsymbol{q}\in\mathbb{R}^{n\times d}$
is linear in $\boldsymbol{q}$ and the divergence operator is defined as
$\nabla\cdot(\mathbf{F}\boldsymbol{q})=\sum_{i=1}^{d}{(A_{i}\boldsymbol{q})}_{x_{i}}$
with matrix functions $A_{i}:\Omega\to\mathbb{R}^{n\times n}$, where the
indices denote partial differentiation with respect to $x_{i}$. Together with
(1a), we assume the boundary and initial conditions $\displaystyle
B\boldsymbol{q}(\boldsymbol{x},t)$
$\displaystyle=\boldsymbol{g}(\boldsymbol{x},t)\qquad$
$\displaystyle\boldsymbol{x}\in\Gamma,\>t\in(0,T),$ (1b)
$\displaystyle\boldsymbol{q}(\boldsymbol{x},0)$
$\displaystyle=\boldsymbol{q}_{0}(\boldsymbol{x})\qquad$
$\displaystyle\boldsymbol{x}\in\Omega.$ (1c)

Here, $B:\Gamma\to\mathbb{R}^{l\times n}$ is a matrix function that takes into
account that boundary conditions can only be prescribed on inflow
characteristics. Under these conditions, (1) has a unique solution
$\boldsymbol{q}$ in a proper space $Q$ GustafssonKreissOliger95 .

We target problems, in which the flux, the right hand side, or the boundary or
initial condition data in (1) depend on parameters $\mathbf{c}$ from a space
$U$. These parameters can either be finite-dimensional, i.e.,
$\mathbf{c}=(c_{1},\ldots,c_{k})$ with $k\geq 1$, or infinite-dimensional,
e.g., a function $\mathbf{c}=\mathbf{c}(\boldsymbol{x})$. Examples for
functions $\mathbf{c}$ are material parameters such as the wave speed, or the
right-hand side forcing in (1a).

Our main interest are inverse and estimation problems, and optimal control
problems governed by hyperbolic systems of the form (1). This leads to
optimization problems of the form

$\min_{\mathbf{c},\boldsymbol{q}}\tilde{\mathcal{J}}(\mathbf{c},\boldsymbol{q})\quad\text{subject
to }~{}\eqref{eq:hyper},$ (2)

where $\tilde{\mathcal{J}}$ is a cost function that depends on the parameters
$\mathbf{c}$ and on the state $\boldsymbol{q}$. The parameters $\mathbf{c}$
may be restricted to an admissible set $U_{a\\!d}\subset U$ for instance to
incorporate bound constraints. If $U_{a\\!d}$ is chosen such that for each
$\mathbf{c}\in U_{a\\!d}$ the state equation (1) admits a unique solution
$\boldsymbol{q}:=\mathcal{S}(\mathbf{c})$ (where $\mathcal{S}$ is the solution
operator for the hyperbolic system), then (2) can be written as an
optimization problem in $\mathbf{c}$ only, namely

$\min_{\mathbf{c}\in
U_{a\\!d}}\mathcal{J}(\mathbf{c}):=\tilde{\mathcal{J}}(\mathbf{c},\mathcal{S}(\mathbf{c})).$
(3)

Existence and (local) uniqueness of solutions to (2) and (3) depend on the
form of the cost function $\tilde{\mathcal{J}}$, properties of the solution
and parameter spaces and of the hyperbolic system and have to be studied on a
case-to-case basis (we refer, for instance, to BorziSchulz12 ; Gunzburger03 ;
Lions85 ; Troltzsch10 ). Our main focus is not the solution of the
optimization problem (3), but the computation of derivatives of $\mathcal{J}$
with respect to $\mathbf{c}$, and the interplay of this derivative computation
with the spatial and temporal discretization of the hyperbolic system (1).
Gradients (and second derivatives) of $\mathcal{J}$ are important to solve (3)
efficiently, and can be used for studying parameter sensitivities or
quantifying the uncertainty in the solution of inverse problems Bui-
ThanhGhattasMartinEtAl13 .

###  2.2 Compatibility of boundary conditions

To ensure the existence of a solution to the adjoint equation, compatibility
conditions between boundary terms in the cost function $\tilde{\mathcal{J}}$,
the boundary operator $B$ in (1b) and the operator $\mathbf{F}$ in (1a) must
hold. We consider cost functions of the form

$\tilde{\mathcal{J}}(\mathbf{c},\boldsymbol{q})=\int_{0}^{T}\\!\\!\\!\int_{\Omega}j_{\Omega}(\boldsymbol{q})\,dx\,dt+\int_{0}^{T}\\!\\!\\!\int_{\Gamma}j_{\Gamma}(C\boldsymbol{q})\,dx\,dt+\int_{\Omega}j_{T}(\boldsymbol{q}(T))\,dx,$
(4)

where $j_{\Omega}:\mathbb{R}^{n}\to\mathbb{R}$,
$j_{\Gamma}:\mathbb{R}^{m}\to\mathbb{R}$ and
$J_{T}:\mathbb{R}^{n}\to\mathbb{R}$ are differentiable, and
$C:\Gamma\to\mathbb{R}^{m\times n}$ is a matrix-valued function. We denote the
derivatives of the functional under the integrals by
$j_{\Omega}^{\prime}(\cdot)$, $j_{\Gamma}^{\prime}(\cdot)$ and
$j_{T}^{\prime}(\cdot)$. The Fréchet derivative of $\mathcal{J}$ with respect
to $\boldsymbol{q}$ in a direction $\tilde{}\boldsymbol{q}$ is given by

$\tilde{\mathcal{J}}_{\boldsymbol{q}}(\mathbf{c},\boldsymbol{q})(\tilde{}\boldsymbol{q})=\int_{0}^{T}\\!\\!\\!\int_{\Omega}j^{\prime}_{\Omega}(\boldsymbol{q})\tilde{}\boldsymbol{q}\,dx\,dt+\int_{0}^{T}\\!\\!\\!\int_{\Gamma}j^{\prime}_{\Gamma}(C\boldsymbol{q})C\tilde{}\boldsymbol{q}\,dx\,dt\\\
+\int_{\Omega}j^{\prime}_{T}(\boldsymbol{q}(T))\tilde{}\boldsymbol{q}(T)\,dx.$
(5)

The boundary operators $B$ and $C$ must be compatible in the sense discussed
next. Denoting the outward pointing normal along the boundary $\Gamma$ by
$\boldsymbol{n}={(n_{1},\ldots,n_{d})}^{T}$, we use the decomposition

$A:=\sum_{i=1}^{d}n_{i}A_{i}=L^{T}\operatorname{diag}(\lambda_{1},\ldots,\lambda_{n})L,$
(6)

with $L\in\mathbb{R}^{n\times n}$ and $\lambda_{1}\geq\ldots\geq\lambda_{n}$.
Note that $L^{-1}=L^{T}$ if $A$ is symmetric. The positive eigenvalues
$\lambda_{1},\ldots,\lambda_{s}$, correspond to the $s\geq 0$ incoming
characteristics, and the negative eigenvalues
$\lambda_{n-m+1},\ldots,\lambda_{n}$ to the $m\geq 0$ outgoing
characteristics. Here, we allow for the zero eigenvalues
$\lambda_{s+1}=\cdots=\lambda_{n-m}=0$. To ensure well-posedness of the
hyperbolic system $\eqref{eq:hyper}$, the initial values of $\boldsymbol{q}$
can only be specified along incoming characteristics. The first $s$ rows
corresponding to incoming characteristics can be identified with the boundary
operator $B$ in (1b). To guarantee well-posedness of the adjoint equation, $C$
has to be chosen such that the cost functional $\tilde{\mathcal{J}}$ only
involves boundary measurements for outgoing characteristics. These correspond
to the rows of $L$ with negative eigenvalues and have to correspond to the
boundary operator $C$. It follows from (6) that

$A=\begin{bmatrix}B\\\ O\\\ C\end{bmatrix}^{-1}\begin{pmatrix}\lambda_{1}\\\
&\ddots\\\ &&\lambda_{n}\end{pmatrix}\begin{bmatrix}B\\\ O\\\
C\end{bmatrix}=\mathrm{bar}C^{T}B-\mathrm{bar}B^{T}C,$ (7)

where $O\in\mathbb{R}^{(n-s-m)\times n}$ and
$\mathrm{bar}B\in\mathbb{R}^{l\times n}$, $\mathrm{bar}C\in\mathbb{R}^{m\times
n}$ are derived properly. If $A$ is symmetric, then
$\mathrm{bar}C^{T}=B^{T}\operatorname{diag}(\lambda_{1},\ldots,\lambda_{S})$,
and
$\mathrm{bar}B^{T}=-C^{T}\operatorname{diag}(\lambda_{n-m+1},\ldots,\lambda_{n})$.
As will be shown in the next section, the matrix $\mathrm{bar}B$ is the
boundary condition matrix for the adjoint equation. For a discussion of
compatibility between the boundary term in a cost functional and hyperbolic
systems in a more general context we refer to GilesPierce97 ; Hartmann07 ;
AlexeSandu10 . In the next section, we formally derive the infinite-
dimensional adjoint system and derivatives of the cost functional
$\tilde{\mathcal{J}}$.

###  2.3 Infinite-dimensional derivatives

For simplicity, we assume that only the flux $\mathbf{F}$ (i.e., the matrices
$A_{1},A_{2},A_{3}$) depend on $\mathbf{c}$, but $B$, $C$, $\boldsymbol{f}$
and $\boldsymbol{q}_{0}$ do not depend on the parameters $\mathbf{c}$. We use
the formal Lagrangian method Troltzsch10 ; BorziSchulz12 , in which we
consider $\mathbf{c}$ and $\boldsymbol{q}$ as independent variables and
introduce the Lagrangian function

$\displaystyle\mathscr{L}(\mathbf{c},\boldsymbol{q},\boldsymbol{p}):=\tilde{\mathcal{J}}(\mathbf{c},\boldsymbol{q})+\int_{0}^{T}\\!\\!\int_{\Omega}{\left(\boldsymbol{q}_{t}+\nabla\cdot(\mathbf{F}\boldsymbol{q})-\boldsymbol{f},\boldsymbol{p}\right)}_{W}\,dx\,dt$
(8)

with $\boldsymbol{p},\boldsymbol{q}\in Q$, where $\boldsymbol{q}$ satisfies
the boundary and initial conditions (1b) and (1c), $\boldsymbol{p}$ satisfies
homogeneous versions of these conditions, and $\mathbf{c}\in U_{a\\!d}$. Here,
${(\cdot\,,\cdot)}_{W}$ denotes a $W$-weighted inner product in
$\mathbb{R}^{n}$, with a symmetric positive definite matrix
$W\in\mathbb{R}^{n\times n}$ (which may depend on $\boldsymbol{x}$). The
matrix $W$ can be used to make a hyperbolic system symmetric with respect to
the $W$-weighted inner product, as for instance in the acoustic and coupled
elastic-acoustic wave examples discussed in Sections 3.2 and 3.4. In
particular, this gives the adjoint equation a form very similar to the state
equation. If boundary conditions depend on the parameters $\mathbf{c}$, they
must be enforced weakly through a Lagrange multiplier in the Lagrangian
function and cannot be added in the definition of the solution space for
$\boldsymbol{q}$. For instance, if the boundary operator $B=B(\mathbf{c})$
depends on $\mathbf{c}$, the boundary condition (1b) must be enforced weakly
through a Lagrange multiplier, amounting to an additional term in the
Lagrangian functional (8).

Following the Lagrangian approach Troltzsch10 ; BorziSchulz12 , the gradient
of $\mathcal{J}$ coincides with the gradient of $\mathscr{L}$ with respect to
$\mathbf{c}$, provided all variations of $\mathscr{L}$ with respect to
$\boldsymbol{q}$ and $\boldsymbol{p}$ vanish. Requiring that variations with
respect to $\boldsymbol{p}$ vanish, we recover the state equation. Variations
with respect to $\boldsymbol{q}$ in directions $\tilde{}\boldsymbol{q}$, that
satisfy homogeneous versions of the initial and boundary conditions (1b) and
(1c), result in

$\displaystyle\mathscr{L_{\boldsymbol{q}}(\mathbf{c},\boldsymbol{q},\boldsymbol{p})}(\tilde{}\boldsymbol{q})$
$\displaystyle=\tilde{\mathcal{J}}_{\boldsymbol{q}}(\mathbf{c},\boldsymbol{q})(\tilde{}\boldsymbol{q})-\int_{0}^{T}\\!\\!\int_{\Omega}\left(\boldsymbol{p}_{t},W\tilde{}\boldsymbol{q}\right)+\left(\mathbf{F}\tilde{}\boldsymbol{q},\nabla(W\boldsymbol{p})\right)\,dx\,dt$
$\displaystyle+\int_{\Omega}\left(\tilde{}\boldsymbol{q}(T),W\boldsymbol{p}(T)\right)+\int_{0}^{T}\\!\\!\int_{\Gamma}\left(\boldsymbol{n}\cdot\mathbf{F}\tilde{}\boldsymbol{q},W\boldsymbol{p}\right)\,dx\,dt,$

where we have used integration by parts in time and space. As will be
discussed in Section 2.5, integration by parts can be problematic when
integrals are approximated using numerical quadrature and should be avoided to
guarantee exact computation of discrete derivatives. In this section, we
assume continuous functions $\boldsymbol{q},\boldsymbol{p}$ and exact
computation of integrals. Since $\boldsymbol{n}\cdot\mathbf{F}=A$, (7) implies
that

$\int_{0}^{T}\\!\\!\int_{\Gamma}(\boldsymbol{n}\cdot\mathbf{F}\tilde{}\boldsymbol{q},W\boldsymbol{p})\,dx\,dt=\int_{0}^{T}\\!\\!\int_{\Gamma}(C\tilde{}\boldsymbol{q},\mathrm{bar}BW\boldsymbol{p})-(B\tilde{}\boldsymbol{q},\mathrm{bar}CW\boldsymbol{p})\,dx\,dt.$
(9)

Using the explicit form of the cost given in (5), and that all variations with
respect to arbitrary $\tilde{}\boldsymbol{q}$ that satisfy
$B\tilde{}\boldsymbol{q}=0$ must vanish, we obtain

$\displaystyle W\boldsymbol{p}_{t}+\mathbf{F}^{\star}\nabla(W\boldsymbol{p})$
$\displaystyle=j^{\prime}_{\Omega}(\boldsymbol{q})\qquad$
$\displaystyle\text{\ on\ }\Omega\times(0,T)$ (10a)
$\displaystyle\mathrm{bar}BW\boldsymbol{p}(\boldsymbol{x},t)$
$\displaystyle=-j^{\prime}_{\Gamma}(C\boldsymbol{q}(\boldsymbol{x},t))\qquad$
$\displaystyle\boldsymbol{x}\in\Gamma,\>t\in(0,T),$ (10b) $\displaystyle
W\boldsymbol{p}(\boldsymbol{x},T)$
$\displaystyle=-j^{\prime}_{T}(\boldsymbol{q}(\boldsymbol{x},T))\qquad$
$\displaystyle\boldsymbol{x}\in\Omega.$ (10c)

Here, $\mathbf{F}^{\star}$ is the adjoint of $\mathbf{F}$ with respect to the
Euclidean inner product. Note that the adjoint system (10) is a final value
problem and thus is usually solved backwards in time. Note that, differently
from the state system, the adjoint system is not in conservative form.

Next we compute variations of $\mathscr{L}$ with respect to the parameters
$\mathbf{c}$ and obtain for variations $\tilde{}\mathbf{c}$ that

$\displaystyle\mathscr{L}_{\mathbf{c}}(\mathbf{c},\boldsymbol{q},\boldsymbol{p})(\tilde{}\mathbf{c})$
$\displaystyle=\int_{0}^{T}\\!\\!\\!\int_{\Omega}{\left(\nabla\cdot(\mathbf{F}_{\mathbf{c}}(\tilde{}\mathbf{c})\boldsymbol{q}),\boldsymbol{p}\right)}_{W}=\int_{0}^{T}\\!\\!\\!\int_{\Omega}\sum_{i=1}^{d}{\left({({A_{i}}_{\mathbf{c}}(\tilde{}\mathbf{c})\boldsymbol{q})}_{x_{i}},\boldsymbol{p}\right)}_{W}$
(11a)

Since $\boldsymbol{q}$ and $\boldsymbol{p}$ are assumed to solve the state and
adjoint system,

$\mathcal{J}(\mathbf{c})(\tilde{}\mathbf{c})=\mathscr{L}_{\mathbf{c}}(\mathbf{c},\boldsymbol{q},\boldsymbol{p})(\tilde{}\mathbf{c}),\
\text{where}\ \boldsymbol{q}\ \text{solves~{}\eqref{eq:hyper}}\ \text{and}\
\boldsymbol{p}\ \text{solves~{}\eqref{eq:adjhyper}}.$ (12)

Next, we present the dG discretization of the hyperbolic system (1) and
discuss the interaction between discretization and the computation of
derivatives.

###  2.4 Discontinuous Galerkin discretization

For the spatial discretization of hyperbolic systems such as (1), the
discontinuous Galerkin (dG) method has proven to be a favorable choice. In the
dG method, we divide the domain $\Omega$ into disjoint elements
${\Omega^{e}}$, and use polynomials to approximate $\boldsymbol{q}$ on each
element ${\Omega^{e}}$. The resulting approximation space is denoted by
$Q^{h}$, and elements $\boldsymbol{q}_{h}\in Q^{h}$ are polynomial on each
element, and discontinuous across elements. Using test functions
$\boldsymbol{p}_{h}\in Q^{h}$, dG discretization in space implies that for
each element ${\Omega^{e}}$

$\begin{split}\int_{{\Omega^{e}}}(\frac{\partial}{\partial
t}\boldsymbol{q}_{h},W\boldsymbol{p}_{h})&-(\mathbf{F}\boldsymbol{q}_{h},\nabla(W\boldsymbol{p}_{h}))\,d\boldsymbol{x}+\\\
&\int_{\Gamma^{e}}\boldsymbol{n}^{-}\cdot({(\mathbf{F}\boldsymbol{q}_{h})}^{\dagger},W\boldsymbol{p}_{h}^{-})\,d\boldsymbol{x}=\int_{{\Omega^{e}}}(\boldsymbol{f},W\boldsymbol{p}_{h})\,d\boldsymbol{x}\end{split}$
(13)

for all times $t\in(0,T)$. Here, $(\cdot\,,\cdot)$ is the inner product in
$\mathbb{R}^{n}$ and $\mathbb{R}^{d\times n}$ and the symmetric and positive
definite matrix $W$ acts as a weighting matrix in this inner product.
Furthermore, ${(\mathbf{F}\boldsymbol{q}_{h})}^{\dagger}$ is the numerical
flux, which connects adjacent elements. The superscript “$-$” denotes that the
inward values are chosen on $\Gamma^{e}$, i.e., the values of the
approximation on $\Omega^{e}$; the superscript “$+$” denotes that the outwards
values are chosen, i.e., the values of an element $\Omega^{e^{\prime}}$ that
is adjacent to $\Omega^{e}$ along the shared boundary $\Gamma^{e}$. Here,
$\boldsymbol{n}^{-}$ is the outward pointing normal on element $\Omega^{e}$.
The formulation (13) is often referred to as the _weak form_ of the dG
discretization HesthavenWarburton08 ; Kopriva09 . The corresponding _strong
form_ dG discretization is obtained by element-wise integration by parts in
space in (13), resulting in

$\begin{split}\int_{{\Omega^{e}}}(\frac{\partial}{\partial
t}\boldsymbol{q}_{h},W\boldsymbol{p}_{h})&+(\nabla\cdot\mathbf{F}\boldsymbol{q}_{h},W\boldsymbol{p}_{h})\,d\boldsymbol{x}-\\\
&\int_{\Gamma^{e}}\boldsymbol{n}^{-}\cdot(\mathbf{F}\boldsymbol{q}_{h}^{-}-{(\mathbf{F}\boldsymbol{q}_{h})}^{\dagger},W\boldsymbol{p}_{h}^{-})\,d\boldsymbol{x}=\int_{{\Omega^{e}}}(\boldsymbol{f},W\boldsymbol{p}_{h})\,d\boldsymbol{x}\end{split}$
(14)

for all $t\in(0,T)$. To find a solution to the optimization problem (3),
derivatives of $\mathcal{J}$ with respect to the parameters $\mathbf{c}$ must
be computed. There are two choices for computing derivatives, namely deriving
expressions for the derivatives of the continuous problem (3), and then
discretizing these equations, or first discretizing the problem and then
computing derivatives of this fully discrete problem. If the latter approach
is taken, i.e., the discrete adjoints are computed, the question arises
weather the discrete adjoint equation is an approximation of the continuous
adjoint and if the discrete adjoint equation is again a dG discretization.
Moreover, what are the consequences of choosing the weak or the strong form
(13) or (14)? A sketch for the different combinations of discretization and
computation of derivatives is also shown in Figure 1.

###  2.5 Influence of numerical quadrature

In a numerical implementation, integrals are often approximated using
numerical quadrature. A fully discrete approach has to take into account the
resulting quadrature error; in particular, integration by parts can incur an
error in combination with numerical quadrature. Below, we first discuss
implications of numerical quadrature in space and then comment on numerical
integration in time. In our example problems in Section 3, we use integral
symbols to denote integration in space and time, but do not assume exact
integration. In particular, we avoid integration by parts or highlight when
integration by parts is used.

${\mathcal{J}}$${\mathcal{J}^{dG}_{h}}$${\mathcal{J}^{dG}_{h,k}}$${\mathcal{J}^{\prime}}$${\mathcal{J}^{dG^{\prime}}_{h}}$${\mathcal{J}^{dG^{\prime}}_{h,k}}$gradgradhgradh,kdG
in spacetime-discretizationdiscretize in spacetime-discretization Figure 1:
Sketch to illustrate the relation between discretization of the problem
(horizontal arrows, upper row), computation of the gradient with respect to
the parameters $\mathbf{c}$ (vertical arrows) and discretization of the
gradient (horizontal arrows, lower row). The problem discretization (upper
row) requires discretization of the state equation and of the cost functional
in space (upper left horizontal arrow) and in time (upper right horizontal
arrow). The vertical arrows represent the Frèchet derivatives of $\mathcal{J}$
(left), of the semidiscrete cost $\mathcal{J}^{dG}_{h}$ (middle) and the fully
discrete cost $\mathcal{J}^{dG}_{h,k}$ (right). The discretization of the
gradient (bottom row) requires space (left arrow) and time (right arrow)
discretization of the state equation, the adjoint equation and the expression
for the gradient. Most of our derivations follow the fully discrete approach,
i.e., the upper row and right arrows; The resulting discrete expressions are
then interpreted as discretizations of the corresponding continuous equations
derived by following the vertical left arrow.

####  2.5.1 Numerical integration in space

The weak form (13) and the strong form (14) of dG are equivalent provided
integrals are computed exactly and, as a consequence, integration by parts
does not result in numerical error. If numerical quadrature is used, these
forms are only numerically equivalent under certain conditions
KoprivaGassner10 ; in general, they are different. To compute fully discrete
gradients, we thus avoid integration by parts in space whenever possible. As a
consequence, if the weak form for the state equation is used, the adjoint
equation is in strong form; this is illustrated and further discussed in
Section 3.

####  2.5.2 Numerical integration in time

We use a method-of-lines approach, that is, from the dG discretization in
space we obtain a continuous-in-time system of ordinary differential equations
(ODEs), which is then discretized by a Runge–Kutta method. Discrete adjoint
equations and the convergence to their continuous counterparts for systems of
ODEs discretized by Runge–Kutta methods have been studied, for instance in
Hager00 ; Walther07 . For the time-discretization we build on these results.

An alternative approach to discretize in time is using a finite element method
for the time discretization, which allows a fully variational formulation of
the problem in space-time; we refer, for instance to BeckerMeidnerVexler07
for this approach applied to parabolic optimization problems.

In both approaches, the computation of derivatives requires the entire time
history of both the state and the adjoint solutions. For realistic application
problems, storing this entire time history is infeasible, and storage
reduction techniques, also known as checkpointing strategies have to be
employed. These methods allow to trade storage against computation time by
storing the state solution only at certain time instances, and then
recomputing it as needed when solving the adjoint equation and computing the
gradient BeckerMeidnerVexler07 ; GriewankWalther00 ; GriewankWalther08 .

##  3 Example problems

The purpose of this section is to illustrate the issues discussed in the
previous sections on example problems. We present examples with increasing
complexity and with parameters entering differently in the hyperbolic systems.
First, in Section 3.1, we derive expressions for derivatives of a functional
that depends on the solution of the one-dimensional advection equation. Since
linear conservation laws can be transformed in systems of advection equations,
we provide extensive details for this example. In particular, we discuss
different numerical fluxes. In Section 3.2, we compute expressions for the
derivatives of a functional with respect to the local wave speed in an
acoustic wave equation. This is followed by examples in which we compute
derivatives with respect to a boundary forcing in Maxwell’s equation (Section
3.3) and derivatives with respect to the primary and secondary wave speeds in
the coupled acoustic-elastic wave equation (Section 3.4). At the end of each
example, we summarize our observations in remarks.

Throughout this section, we use the dG discretization introduced in Section
2.4 and denote the finite dimensional dG solution spaces by $P^{h}$ and
$Q^{h}$. These spaces do not include the boundary conditions which are usually
imposed weakly through the numerical flux in the dG method, and they also do
not include initial/final time conditions, which we specify explicitly.
Functions in these dG spaces are smooth (for instance polynomials) on each
element $\Omega^{e}$ and discontinuous across the element boundaries
$\partial{\Omega^{e}}$. As before, for each element we denote the inward value
with a superscript “-” and the outward value with superscript “+”. We use the
index $h$ to denote discretized fields and denote by $[\\![\cdot]\\!]$ the
jump, by $[\cdot]$ the difference and by $\\{\\!\\!\\{\cdot\\}\\!\\!\\}$ the
mean value at an element interface $\partial{\Omega^{e}}$. To be precise, for
a scalar dG-function $u_{h}$, these are defined as
$[\\![u_{h}]\\!]=u_{h}^{-}\boldsymbol{n}^{-}+u_{h}^{+}\boldsymbol{n}^{+}$,
$[u_{h}]=u_{h}^{-}-u_{h}^{+}$ and
$\\{\\!\\!\\{u_{h}\\}\\!\\!\\}=(u_{h}^{-}+u_{h}^{+})/2$. Likewise, for a
vector $\boldsymbol{v}_{h}$ we define
$[\\![\boldsymbol{v}_{h}]\\!]=\boldsymbol{v}_{h}^{-}\cdot\boldsymbol{n}^{-}+\boldsymbol{v}_{h}^{+}\cdot\boldsymbol{n}^{+}$,
$[\boldsymbol{v}_{h}]=\boldsymbol{v}_{h}^{-}-\boldsymbol{v}_{h}^{+}$ and
$\\{\\!\\!\\{\boldsymbol{v}_{h}\\}\\!\\!\\}=(\boldsymbol{v}_{h}^{-}+\boldsymbol{v}_{h}^{+})/2$
and for a second-order tensor $\SS_{h}$ we have
$[\\![\SS_{h}]\\!]=\SS_{h}^{-}\boldsymbol{n}^{-}+\SS_{h}^{+}\boldsymbol{n}^{+}$.
The domain boundary $\partial\Omega$ is denoted by $\Gamma$.

Throughout this section, we use the usual symbol to denote integrals, but we
do not assume exact quadrature. Rather, integration can be replaced by a
numerical quadrature rule and, as a consequence, the integration by parts
formula does not hold exactly. To avoid numerical errors when using numerical
quadrature, we thus avoid integration by parts in space or point out when
integration by parts is used. Since our focus is on the spatial dG
discretization, we do not discretize the problem in time and assume exact
integration in time.

###  3.1 One-dimensional advection equation

We consider the one-dimensional advection equation on the spatial domain
$\Omega=(x_{l},x_{r})\subset\mathbb{R}$. We assume a spatially varying,
continuous positive advection velocity $a(x)\geq a_{0}>0$ for $x\in\Omega$ and
a forcing $f(x,t)$ for $(x,t)\in\Omega\times(0,T)$. The advection equation
written in conservative form is given by

$\displaystyle u_{t}+{(au)}_{x}$ $\displaystyle=f\quad$ $\displaystyle\text{\
on\ }\Omega\times(0,T),$ (15a) with the initial condition $\displaystyle
u(x,0)$ $\displaystyle=u_{0}(x)\quad$ $\displaystyle\text{\ for\ }x\in\Omega,$
(15b) and, since $a>0$, the inflow boundary is $\Gamma_{l}:=\\{x_{l}\\}$,
where we assume $\displaystyle u(x_{l},t)$ $\displaystyle=u_{l}(t)\quad$
$\displaystyle\text{\ for\ }t\in(0,T).$ (15c)

The discontinuous Galerkin (dG) method for the numerical solution of (15a) in
strong form is: Find $u_{h}\in P^{h}$ with $u_{h}(x,0)=u_{0}(x),x\in\Omega$
such that for all test functions $p_{h}\in P^{h}$ holds

$\int_{\Omega}(u_{h,t}+{(au_{h})}_{x}-f)p_{h}\,dx=\sum_{e}\int_{\partial{\Omega^{e}}}n^{-}\left(au_{h}^{-}-{(au_{h})}^{\dagger}\right)p_{h}^{-}\,dx$
(16)

for all $t\in(0,T)$. Here, $a\in U$ can be an infinite-dimensional continuous
function, or a finite element function. For $\alpha\in[0,1]$, the numerical
flux on the boundary is replaced by the numerical flux,
${(au_{h})}^{\dagger}$, given by

${(au_{h})}^{\dagger}=a\\{\\!\\!\\{u_{h}\\}\\!\\!\\}+\frac{1}{2}\left\lvert
a\right\rvert(1-\alpha)[\\![u_{h}]\\!].$ (17)

This is a central flux for $\alpha=1$, and an upwind flux for $\alpha=0$. Note
that in one spatial dimension, the outward normal $n$ is $-1$ and $+1$ on the
left and right side of ${\Omega^{e}}$, respectively. Since $a$ is assumed to
be continuous on $\Omega$, we have $a^{-}=a^{+}:=a$. Moreover, since $a$ is
positive, we neglect the absolute value in the following.

As is standard practice HesthavenWarburton08 ; Kopriva09  we incorporate the
boundary conditions weakly through the numerical flux by choosing the
“outside” values $u_{h}^{+}$ as $u_{h}^{+}=u_{l}$ for $x=x_{l}$ and
$u_{h}^{+}=u_{h}^{-}$ for $x=x_{r}$ for the computation of
${(au_{h})}^{\dagger}$. This implies that
$[\\![u_{h}]\\!]=n^{-}(u_{h}^{-}-u_{l})$ on $\Gamma_{l}$ and
$[\\![u_{h}]\\!]=0$ on the outflow boundary $\Gamma_{r}:=\\{x_{r}\\}$. For
completeness, we also provide the dG discretization of (15a) in weak form:
Find $u_{h}\in P^{h}$ with $u_{h}(x,0)=u_{0}(x),x\in\Omega$ such that for all
$p_{h}\in P^{h}$ holds

$\int_{\Omega}(u_{h,t}-f)p_{h}-au_{h}p_{h,x}\,dx=-\sum_{e}\int_{\partial{\Omega^{e}}}n^{-}{(au_{h})}^{\dagger}p_{h}^{-}\,dx$

for all $t\in(0,T)$, which is found by element-wise integration by parts in
(16).

Adding the boundary contributions in (16) from adjacent elements
${\Omega^{e}}$ and ${\Omega^{e^{\prime}}}$ to their shared point
$\partial{\Omega^{e}}\cap\partial{\Omega^{e^{\prime}}}$, we obtain

$\frac{1}{2}\left\lvert
a\right\rvert(\alpha-1)[\\![u_{h}]\\!][\\![p_{h}]\\!]+a[\\![u_{h}]\\!]\\{\\!\\!\\{p_{h}\\}\\!\\!\\}.$

Thus, (16) can also be written as

$\int_{\Omega}(u_{h,t}+{(au_{h})}_{x}-f)p_{h}\,dx=\sum_{e}\int_{\partial{\Omega^{e}}\setminus\Gamma}n^{-}\left(a\\{\\!\\!\\{p_{h}\\}\\!\\!\\}+{a}\frac{1}{2}(\alpha-1)[\\![p_{h}]\\!]\right)u_{h}^{-}\,dx\\\
+\int_{\Gamma_{l}}n^{-}\left(a\frac{1}{2}p_{h}^{-}+{a}\frac{1}{2}(\alpha-1)n^{-}p_{h}^{-}\right)(u_{h}^{-}-u_{l})\,dx.$
(18)

We consider an objective functional for the advection velocity $a\in U$ given
by

$\mathcal{J}(a):=\tilde{\mathcal{J}}(a,u_{h})=\int_{0}^{T}\\!\\!\int_{\Omega}j_{\Omega}(u_{h})\,dx\,dt+\int_{0}^{T}\\!\\!\int_{\Gamma}j_{\Gamma}(u_{h})\,dx\,dt+\int_{\Omega}r(a)\,dx,$
(19)

with differentiable functions $j_{\Omega}:\Omega\to\mathbb{R}$,
$j_{\Gamma}:\Gamma\to\mathbb{R}$ and $r:\Omega\to\mathbb{R}$. For illustration
purposes, we define the boundary term in the cost on both, the inflow and the
outflow part of the boundary, and comment on the consequences in Remark 1. To
derive the discrete gradient of $\mathcal{J}$, we use the Lagrangian function
$\mathcal{L}:U\times P^{h}\times P^{h}\rightarrow\mathbb{R}$, which combines
the cost (19) with the dG discretization (18):

$\displaystyle\mathcal{L}(a,u_{h},p_{h}):=$
$\displaystyle\tilde{\mathcal{J}}(a,u_{h})+\int_{0}^{T}\\!\\!\int_{\Omega}(u_{h,t}+{(au_{h})}_{x}-f)p_{h}\,dx\,dt$
$\displaystyle-\sum_{e}\int_{0}^{T}\\!\\!\int_{\partial{\Omega^{e}}\setminus\Gamma}n^{-}\left(a\\{\\!\\!\\{p_{h}\\}\\!\\!\\}+{a}\frac{1}{2}(\alpha-1)[\\![p_{h}]\\!]\right)u_{h}^{-}\,dx\,dt$
$\displaystyle-\int_{0}^{T}\int_{\Gamma_{l}}n^{-}\left(a\frac{1}{2}p_{h}^{-}+{a}\frac{1}{2}(\alpha-1)n^{-}p_{h}^{-}\right)(u_{h}^{-}-u_{l})\,dx\,dt.$

By requiring that all variations with respect to $p_{h}$ vanish, we recover
the state equation (18). Variations with respect to $u_{h}$ in a direction
$\tilde{u}_{h}$, which satisfies the homogeneous initial conditions
$\tilde{u}_{h}(x,0)=0$ for $x\in\Omega$ result in

$\displaystyle\\!\\!\\!\\!\mathcal{L}_{u_{h}}(a,u_{h},p_{h})(\tilde{u}_{h})$
$\displaystyle=$
$\displaystyle\int_{0}^{T}\\!\\!\int_{\Omega}j^{\prime}_{\Omega}(u_{h})\tilde{u}_{h}-p_{h,t}\tilde{u}_{h}+{(a\tilde{u}_{h})}_{x}p_{h}\,dx\,dt+\int_{\Omega}p_{h}(x,T)\tilde{u}_{h}(x,T)\,dx$
$\displaystyle-\sum_{e}\int_{0}^{T}\\!\\!\int_{\partial{\Omega^{e}}\setminus\Gamma}\\!\\!\\!n^{-}\left(a\\{\\!\\!\\{p_{h}\\}\\!\\!\\}+{a}\frac{1}{2}(\alpha-1)[\\![p_{h}]\\!]\right)\tilde{u}_{h}^{-}\,dx\,dt+\int_{0}^{T}\\!\\!\int_{\Gamma}j^{\prime}_{\Gamma}(u_{h})\tilde{u}_{h}\,dx\,dt$
$\displaystyle-\int_{0}^{T}\\!\\!\int_{\Gamma_{l}}n^{-}\left(a\frac{1}{2}p_{h}^{-}+{a}\frac{1}{2}(\alpha-1)n^{-}p_{h}^{-}\right)\tilde{u}_{h}^{-}\,dx\,dt$

Since we require that arbitrary variations with respect to $u_{h}$ must
vanish, $p_{h}$ has to satisfy $p_{h}(x,T)=0$ for all $x\in\Omega$, and

$\int_{\Omega}-p_{h,t}\tilde{u}_{h}+{(a\tilde{u}_{h})}_{x}p_{h}+j_{\Omega}^{\prime}(u_{h})\tilde{u}_{h}\,dx=\sum_{e}\int_{\partial{\Omega^{e}}\setminus\Gamma_{l}}n^{-}{(ap_{h})}^{\dagger}\tilde{u}_{h}^{-}\,dx\\\
+\int_{\Gamma_{l}}n^{-}\left(\frac{a(2-\alpha)}{2}p_{h}^{-}+j_{\Gamma}^{\prime}(u_{h})\right)\tilde{u}_{h}^{-}\,dx$
(20a) for all $t\in(0,T)$ and for all $\tilde{u}_{h}$, with the adjoint flux
${(ap_{h})}^{\dagger}:=a\\{\\!\\!\\{p_{h}\\}\\!\\!\\}+\frac{1}{2}{a}(\alpha-1)[\\![p_{h}]\\!]$
(20b) and with
$p_{h}^{+}:=-\frac{j_{\Gamma}^{\prime}(u_{h})}{a(1-\frac{\alpha}{2})}-\frac{\alpha}{2-\alpha}p_{h}^{-}\quad\text{on}\
\Gamma_{r}.$ (20c)

The outside value $p_{h}^{+}$ in (20c) is computed such that
$n^{-}{(ap_{h})}^{\dagger}=j^{\prime}_{\Gamma}(u_{h})$ on $\Gamma_{r}$. The
equations (20) are the weak form of a discontinuous Galerkin discretization of
the adjoint equation, with flux given by (20b). An element-wise integration by
parts in space in (20a), results in the corresponding strong form of the
discrete adjoint equation:

$\begin{split}\int_{\Omega}\big{(}-p_{h,t}-ap_{h,x}+&j_{\Omega}^{\prime}(u_{h})\big{)}\tilde{u}_{h}\,dx=-\sum_{e}\int_{\partial{\Omega^{e}}\setminus\Gamma_{l}}n^{-}\left(ap_{h}^{-}-{(ap_{h})}^{\dagger}\right)\tilde{u}_{h}^{-}\,dx\\\
&+\int_{\Gamma_{l}}n^{-}\left(-\frac{a\alpha}{2}p_{h}^{-}+j_{\Gamma}^{\prime}(u_{h})\right)\tilde{u}_{h}^{-}\,dx\end{split}$
(21)

Note that this integration by parts is not exact if numerical quadrature is
used. It can be avoided if the dG weak form of the adjoint equation (20) is
implemented directly.

Provided $u_{h}$ and $p_{h}$ are solutions to the state and adjoint equations,
respectively, the gradient of $\mathcal{J}$ with respect to $a$ is found by
taking variations of the Lagrangian with respect to $a$ in a direction
$\tilde{a}$:

$\displaystyle\mathcal{L}_{a}(a,u_{h},p_{h})(\tilde{a})=$
$\displaystyle\int_{\Omega}r^{\prime}(a)\tilde{a}\,dx+\int_{0}^{T}\\!\\!\int_{\Omega}{(\tilde{a}u_{h})}_{x}p_{h}\,dx\,dt$
$\displaystyle-\sum_{e}\int_{0}^{T}\\!\\!\int_{\partial{\Omega^{e}}\setminus\Gamma}n^{-}\tilde{a}\left(\frac{1}{2}(\alpha-1)[\\![p_{h}]\\!]+\\{\\!\\!\\{p_{h}\\}\\!\\!\\}\right)u_{h}^{-}\,dx\,dt$
$\displaystyle-\int_{0}^{T}\\!\\!\int_{\Gamma_{l}}n^{-}\tilde{a}\left(\frac{1}{2}(\alpha-1)n^{-}p_{h}^{-}+\frac{1}{2}p_{h}^{-}\right)(u_{h}^{-}-u_{l})\,dx\,dt.$

Thus, the gradient of $\mathcal{J}$ is given by

${\mathcal{J}}^{\prime}(a)(\tilde{a})=\int_{\Omega}r^{\prime}(a)\tilde{a}\,dx+\int_{0}^{T}\\!\\!\int_{\Omega}G\tilde{a}\,dx\,dt+\sum_{\mathsf{f}}\int_{0}^{T}\\!\\!\int_{\mathsf{f}}g\tilde{a}\,dx\,dt,$
(22)

where $G$ is defined for each element ${\Omega^{e}}$, and $g$ for each inter-
element face as follows:

$\displaystyle G\tilde{a}=p_{h}{(\tilde{a}u_{h})}_{x}\quad\text{and}\quad
g=\begin{cases}\frac{1}{2}(\alpha-1)[\\![u_{h}]\\!][\\![p_{h}]\\!]+[\\![u_{h}]\\!]\\{\\!\\!\\{p_{h}\\}\\!\\!\\}&\
\text{if}\ \mathsf{f}\not\subset\Gamma,\\\ \frac{1}{2}\alpha
p_{h}^{-}(u_{h}^{-}-u_{l})&\ \text{if}\ \mathsf{f}\subset\Gamma_{l},\\\ 0&\
\text{if}\ \mathsf{f}\subset\Gamma_{r}.\end{cases}$

Note that the element boundary jump terms arising in $\mathcal{J}^{\prime}(a)$
are a consequence of using the dG method to discretize the state equation.
While in general these terms do not vanish, they become small as the
discretization resolves the continuous state and adjoint variables. However,
these terms must be taken into account to compute discretely exact gradients.
We continue with a series of remarks:

###### Remark 1

The boundary conditions for the adjoint variable $p_{h}$ that are weakly
imposed through the adjoint numerical flux are the Dirichlet condition
$ap_{h}(x_{r})=-j^{\prime}_{\Gamma}(u)$ at $\Gamma_{r}$, which, due to the
sign change for the advection term in (20) and (21), is an inflow boundary for
the adjoint equation. At the (adjoint) outflow boundary $\Gamma_{l}$, the
adjoint scheme can only be stable if ${j_{\Gamma}}|_{\Gamma_{l}}\equiv 0$.
This corresponds to the discussion from Section 2.2 on the compatibility of
boundary operators. The discrete adjoint scheme is consistent (in the sense
that the continuous adjoint variable $p$ satisfies the discrete adjoint
equation (20)) when $\alpha=0$, i.e., for a dG scheme based on an upwind
numerical flux. Thus, only dG discretizations based on upwind fluxes at the
boundary can be used in adjoint calculus. Hence, in the following we restrict
ourselves to upwind fluxes.

###### Remark 2

While the dG discretization of the state equation is in conservative form, the
discrete adjoint equation is not. Moreover, using dG method in strong form for
the state system, the adjoint system is naturally dG method in weak form (see
(20)), and element-wise integration by parts is necessary to find the adjoint
in strong form (21). Vice versa, using the weak form of dG for the state
equation, the adjoint equation is naturally in strong form. These two forms
can be numerically different if the integrals are approximated through a
quadrature rule for which integration by parts does not hold exactly. In this
case, integration by parts should be avoided to obtain exact discrete
derivatives.

###### Remark 3

The numerical fluxes (17) and (20b) differ by the sign in the upwinding term
only. Thus, an upwinding flux for the state equation becomes a downwinding
flux for the adjoint equation. This is natural since the advection velocity
for the adjoint equation is $-a$, which makes the adjoint numerical flux an
upwind flux for the adjoint equation.

###  3.2 Acoustic wave equation

Next, we derive expressions for the discrete gradient with respect to the
local wave speed in the acoustic wave equation. This is important, for
instance, in seismic inversion using full wave forms Fichtner11 ;
LekicRomanowicz11 ; EpanomeritakisAkccelikGhattasEtAl08 ;
PeterKomatitschLuoEtAl11 . If the dG method is used to discretize the wave
equation (as, e.g., in Bui-ThanhBursteddeGhattasEtAl12 ; Bui-
ThanhGhattasMartinEtAl13 ; CollisObervanBloemenWaanders10 ), the question on
the proper discretization of the adjoint equation and of the expressions for
the derivatives arises. Note that in Section 3.4 we present the discrete
derivatives with respect to the (possibly discontinuous) primary and secondary
wave speeds in the coupled acoustic-elastic wave equation, generalizing the
results presented in this section. However, for better readability we choose
to present this simpler case first and then present the results for the
coupled acoustic-elastic equation in compact form in Section 3.4.

We consider the acoustic wave equation written as first-order system as
follows:

$\displaystyle e_{t}-\nabla\cdot\boldsymbol{v}$ $\displaystyle=0\quad$
$\displaystyle\text{\ on\ }\Omega\times(0,T),$ (23a)
$\displaystyle\rho\boldsymbol{v}_{t}-\nabla(\lambda e)$
$\displaystyle=\boldsymbol{f}\quad$ $\displaystyle\text{\ on\
}\Omega\times(0,T),$ (23b) where $\boldsymbol{v}$ is the velocity, $e$ the
dilatation (trace of the strain tensor), $\rho=\rho(\boldsymbol{x})$ is the
mass density, and $\lambda=c^{2}\rho$, where $c(\boldsymbol{x})$ denotes the
wave speed. Together with (23a) and (23b), we assume the initial conditions
$\displaystyle
e(\boldsymbol{x},0)=e_{0}(\boldsymbol{x}),\>\boldsymbol{v}(\boldsymbol{x},0)=\boldsymbol{v}_{0}(\boldsymbol{x})\quad\text{\
for\ }\boldsymbol{x}\in\Omega,$ (23c) and the boundary conditions
$\displaystyle
e(\boldsymbol{x},t)=e_{\text{bc}}(\boldsymbol{x},t),\>\>\boldsymbol{v}(\boldsymbol{x},t)=\boldsymbol{v}_{\text{bc}}(\boldsymbol{x},t)\quad\text{\
for\ }(\boldsymbol{x},t)\in\Gamma\times(0,T).$ (23d)

Note that the dG method discussed below uses an upwind numerical flux, and
thus the boundary conditions (23) are automatically only imposed at inflow
boundaries. Through proper choice of $e_{\text{bc}}$ and
$\boldsymbol{v}_{\text{bc}}$ classical wave equation boundary conditions can
be imposed, e.g., WilcoxStadlerBursteddeEtAl10 ; FengTengChen07 .

The choice of the dilatation $e$ together with the velocity $\boldsymbol{v}$
in the first order system formulation is motivated from the strain-velocity
formulation used for the coupled elastic and acoustic wave equation in Section
3.4. To write (23a) and (23b) in second-order form, we define the pressure as
$p=-\lambda e$ and obtain the pressure-velocity form as

$\displaystyle p_{t}+\lambda\nabla\cdot\boldsymbol{v}$ $\displaystyle=0\quad$
$\displaystyle\text{\ on\ }\Omega\times(0,T),$
$\displaystyle\rho\boldsymbol{v}_{t}+\nabla p$
$\displaystyle=\boldsymbol{f}\quad$ $\displaystyle\text{\ on\
}\Omega\times(0,T),$

which is equivalent to the second-order formulation

$p_{tt}=\lambda\nabla\cdot\left(\frac{1}{\rho}\nabla
p\right)-\lambda\nabla\cdot\left(\frac{1}{\rho}\boldsymbol{f}\right)\quad\text{\
on\ }\Omega\times(0,T).$ (25)

The strong form dG discretization of (23) is: Find
$(e_{h},\boldsymbol{v}_{h})\in P^{h}\times Q^{h}$ satisfying the initial
conditions (23c) such that for all test functions
$(h_{h},\boldsymbol{w}_{h})\in P^{h}\times Q^{h}$ and for all $t\in(0,T)$
holds:

$\begin{split}&\int_{\Omega}(e_{t,h}-\nabla\cdot\boldsymbol{v}_{h})\lambda
h_{h}\,d\boldsymbol{x}+\int_{\Omega}(\rho\boldsymbol{v}_{t,h}-\nabla(\lambda
e_{h})-\boldsymbol{f})\cdot\boldsymbol{w}_{h}\,d\boldsymbol{x}\\\
&=-\sum_{e}\int_{\partial{\Omega^{e}}}\boldsymbol{n}^{-}\cdot\left(\boldsymbol{v}_{h}^{-}-\boldsymbol{v}_{h}^{\dagger}\right)\lambda
h_{h}^{-}+\left({(\lambda e_{h})}^{-}-{(\lambda
e_{h})}^{\dagger}\right)\boldsymbol{n}^{-}\cdot\boldsymbol{w}_{h}^{-}\,d\boldsymbol{x}.\end{split}$
(26)

Note that above, the inner product used for (23a) is weighted by $\lambda$,
which makes the first-order form of the wave equation a symmetric operator and
also allows for a natural interpretation of the adjoint variables, as shown
below. Assuming that $c$ and $\rho$ are continuous, we obtain the upwind
numerical fluxes:

$\displaystyle\boldsymbol{n}^{-}\cdot\boldsymbol{v}_{h}^{\dagger}$
$\displaystyle=\boldsymbol{n}^{-}\cdot\\{\\!\\!\\{\boldsymbol{v}_{h}\\}\\!\\!\\}-\frac{c}{2}[e_{h}],$
(27a) $\displaystyle{(\lambda e_{h})}^{\dagger}$
$\displaystyle=\lambda\\{\\!\\!\\{e_{h}\\}\\!\\!\\}-\frac{\rho
c}{2}[\\![\boldsymbol{v}_{h}]\\!].$ (27b)

Adding the boundary contributions from two adjacent elements ${\Omega^{e}}$
and ${\Omega^{e^{\prime}}}$ in (26) to a shared edge (in 2D) or face (in 3D),
one obtains

$\lambda[\\![\boldsymbol{v}_{h}]\\!]\\{\\!\\!\\{h_{h}\\}\\!\\!\\}+\frac{c\lambda}{2}[e_{h}][h_{h}]+\frac{\rho
c}{2}[\\![\boldsymbol{v}_{h}]\\!][\\![\boldsymbol{w}_{h}]\\!]+\lambda[\\![e_{h}]\\!]\cdot\\{\\!\\!\\{\boldsymbol{w}_{h}\\}\\!\\!\\}.$
(28)

We compute the discrete gradient with respect to the wave speed $c$ for a cost
functional of the form

$\tilde{\mathcal{J}}(c,\boldsymbol{v}_{h})=\int_{0}^{T}\\!\\!\int_{\Omega}j_{\Omega}(\boldsymbol{v}_{h})\,d\boldsymbol{x}\,dt+\int_{0}^{T}\\!\\!\int_{\Gamma}j_{\Gamma}(\boldsymbol{v}_{h})\,d\boldsymbol{x}\,dt+\int_{\Omega}r(c)\,d\boldsymbol{x},$
(29)

with differentiable functions $j_{\Omega}:\Omega\to\mathbb{R}$,
$j_{\Gamma}:\Gamma\to\mathbb{R}$ and $r:\Omega\to\mathbb{R}$. To ensure
compatibility as discussed in Section 2.2, the boundary term $j_{\Gamma}$ in
(29) can only involve outgoing characteristics. We introduce the Lagrangian
function, use that all its variations with respect to $\boldsymbol{v}$ and $e$
must vanish, and integrate by parts in time $t$, resulting in the following
adjoint equation: Find $(h_{h},\boldsymbol{w}_{h})\in P^{h}\times Q^{h}$
satisfying the finial time conditions $h(\boldsymbol{x},T)=0$,
$\boldsymbol{w}(\boldsymbol{x},T)=0$ for $\boldsymbol{x}\in\Omega$, such that
for all test functions $(\tilde{e}_{h},\tilde{}\boldsymbol{v}_{h})\in
P^{h}\times Q^{h}$ and for all $t\in(0,T)$ holds:

$\begin{split}&\int_{\Omega}-\tilde{e}_{h}\lambda
h_{h,t}-\nabla\cdot\tilde{}\boldsymbol{v}_{h}\lambda
h_{h}-\rho\tilde{}\boldsymbol{v}_{h}\cdot\boldsymbol{w}_{h,t}-\nabla(\lambda\tilde{e}_{h})\cdot\boldsymbol{w}_{h}+j_{\Omega}^{\prime}(\boldsymbol{v}_{h})\cdot\tilde{}\boldsymbol{v}_{h}\,d\boldsymbol{x}\\\
&=-\sum_{e}\int_{\partial{\Omega^{e}}}\left(\boldsymbol{n}^{-}\cdot\boldsymbol{w}_{h}^{\dagger}\right)\lambda\tilde{e}_{h}^{-}+{(\lambda
h_{h})}^{\dagger}\boldsymbol{n}^{-}\cdot\tilde{}\boldsymbol{v}_{h}^{-}\,d\boldsymbol{x}-\int_{\Gamma}j_{\Gamma}^{\prime}(\boldsymbol{v}_{h})\cdot\tilde{}\boldsymbol{v}_{h}\,d\boldsymbol{x},\end{split}$
(30)

where the adjoint numerical fluxes are given by

$\displaystyle\boldsymbol{n}^{-}\cdot\boldsymbol{w}_{h}^{\dagger}$
$\displaystyle=\boldsymbol{n}^{-}\cdot\\{\\!\\!\\{\boldsymbol{w}_{h}\\}\\!\\!\\}+\frac{c}{2}[h_{h}],$
(31a) $\displaystyle{(\lambda h_{h})}^{\dagger}$
$\displaystyle=\lambda\\{\\!\\!\\{h_{h}\\}\\!\\!\\}+\frac{\rho
c}{2}[\\![\boldsymbol{w}_{h}]\\!].$ (31b)

Note that (30) is the weak form of the dG discretization for an acoustic wave
equation, solved backwards in time. This is a consequence of the symmetry of
the differential operator in the acoustic wave equation, when considered in
the appropriate inner product. Comparing (31) and (27) shows that the adjoint
numerical flux (31) is the downwind flux in the adjoint variables for the
adjoint wave equation. The strong dG form corresponding to (30) can be
obtained by element-wise integration in parts in space.

Finally, we present expressions for the derivative of $\mathcal{J}$ with
respect to the wave speed $c$, which are found as variations of the Lagrangian
with respect to $c$. This results in

${\mathcal{J}}^{\prime}(c)(\tilde{c})=\int_{\Omega}r^{\prime}(c)\tilde{c}\,d\boldsymbol{x}+\int_{0}^{T}\\!\\!\int_{\Omega}G\tilde{c}\,d\boldsymbol{x}\,dt+\sum_{\mathsf{f}}\int_{0}^{T}\\!\\!\int_{\mathsf{f}}g\tilde{c}\,d\boldsymbol{x}\,dt,$
(32)

where $G$ is defined on each element ${\Omega^{e}}$, and $g$ for each inter-
element face $\mathsf{f}$ as follows:

$\displaystyle G\tilde{c}$ $\displaystyle=-2\nabla(\rho
c\tilde{c}e_{h})\cdot\boldsymbol{w}_{h}+2\rho
c\tilde{c}(e_{h,t}-\nabla\cdot\boldsymbol{v}_{h})h_{h},$ (33a) $\displaystyle
g$ $\displaystyle=2\rho
c[\\![\boldsymbol{v}_{h}]\\!]\\{\\!\\!\\{h_{h}\\}\\!\\!\\}+\frac{3}{2}\rho
c^{2}[e_{h}][h_{h}]+\frac{1}{2}\rho[\\![\boldsymbol{v}_{h}]\\!][\\![\boldsymbol{w}_{h}]\\!]+2\rho
c[\\![e_{h}]\\!]\cdot\\{\\!\\!\\{\boldsymbol{w}_{h}\\}\\!\\!\\}$ (33b)

Above, $(\boldsymbol{v}_{h},e_{h})$ is the solution of the state equation (26)
and $(\boldsymbol{w}_{h},h_{h})$ the solution of the adjoint equation (30).
Since the state equation (23) is satisfied in the dG sense, (33) simplifies
provided $\rho c\tilde{c}h_{h}\in P^{h}$, or if a quadrature method is used in
which the values of $\rho c\tilde{c}h_{h}$ at the quadrature points coincide
with the values of a function in $P^{h}$ at these points. The latter is, for
instance, always the case when the same nodes are used for the quadrature and
the nodal basis. Then, (33) simplifies to

$\displaystyle G\tilde{c}$ $\displaystyle=-2\nabla(\rho
c\tilde{c}e_{h})\cdot\boldsymbol{w}_{h},$ (34a) $\displaystyle g$
$\displaystyle=\frac{1}{2}\rho
c^{2}[e_{h}][h_{h}]+\frac{1}{2}\rho[\\![\boldsymbol{v}_{h}]\\!][\\![\boldsymbol{w}_{h}]\\!]+2\rho
c[\\![e_{h}]\\!]\cdot\\{\\!\\!\\{\boldsymbol{w}_{h}\\}\\!\\!\\}.$ (34b)

As for the one-dimensional advection problem (see Remark 3), the upwind flux
in the state equation becomes a downwind flux in the adjoint equation, and
thus an upwind flux for the adjoint equation when solved backwards in time.

###### Remark 4

As in the advection example, the discrete gradient has boundary contributions
that involve jumps of the dG variables at the element boundaries (see (33b)
and (34b)). These jumps are at the order of the dG approximation error and
thus tend to zero as the dG solution converges to the continuous solution
either through mesh refinement or improvement of the approximation on each
element.

###  3.3 Maxwell’s equations

Here we derive expressions for the discrete gradient with respect to the
current density in Maxwell’s equations (specifically boundary current density
in our case). This can be used, for instance, in the determination and
reconstruction of antennas from boundary field measurements Nicaise00  and
controlling electromagnetic fields using currents Lagnese89 ; Yousept12 . The
time-dependent Maxwell’s equations in a homogeneous isotropic dielectric
domain $\Omega\subset\mathbb{R}^{3}$ is given by:

$\displaystyle\mu\boldsymbol{H}_{t}$
$\displaystyle=-\nabla\times\boldsymbol{E}$ $\displaystyle\qquad\text{on}\
\Omega\times(0,T),$ (35a) $\displaystyle\epsilon\boldsymbol{E}_{t}$
$\displaystyle=\phantom{-}\nabla\times\boldsymbol{H}$
$\displaystyle\qquad\text{on}\ \Omega\times(0,T),$ (35b)
$\displaystyle\nabla\cdot\boldsymbol{H}$ $\displaystyle=0$
$\displaystyle\qquad\text{on}\ \Omega\times(0,T),$ (35c)
$\displaystyle\nabla\cdot\boldsymbol{E}$ $\displaystyle=0$
$\displaystyle\qquad\text{on}\ \Omega\times(0,T),$ (35d) where
$\boldsymbol{E}$ is the electric field and $\boldsymbol{H}$ is the magnetic
field. Moreover, $\mu$ is the permeability and $\epsilon$ is the permittivity,
which can both be discontinuous across element interfaces. The impedance $Z$
and the conductance $Y$ of the material are defined as
$Z=\frac{1}{Y}=\sqrt{\frac{\mu}{\epsilon}}$. Note that we follow a standard
notation for Maxwell’s equation, in which the vectors $\boldsymbol{H}$ and
$\boldsymbol{E}$ are denoted by bold capital letters. Together with equations
(35a)–(35d), we assume the initial conditions
$\displaystyle\boldsymbol{E}(\boldsymbol{x},0)=\boldsymbol{E}_{0}(\boldsymbol{x}),\>\boldsymbol{H}(\boldsymbol{x},0)=\boldsymbol{H}_{0}(\boldsymbol{x})\quad\text{on}\
\Omega,$ (35e) and boundary conditions
$\displaystyle\boldsymbol{n}\times\boldsymbol{H}=-\boldsymbol{J}_{s}\quad\text{on}\
\Gamma.$ (35f) This classic boundary condition can be converted to equivalent
inflow characteristic boundary conditions TengLinChangEtAl08 . Here,
$\boldsymbol{J}_{s}(\boldsymbol{x},t)$ is a spatially (and possibly time-
dependent) current density flowing tangentially to the boundary.

If the initial conditions satisfy the divergence conditions (35c) and (35d),
the time evolved solution will as well HesthavenWarburton02 . Thus, the
divergence conditions can be regarded as a consistency condition on the
initial conditions. We consider a dG discretization of Maxwell’s equations hat
only involves equations (35a) and (35b) explicitly. The dG solution then
satisfies the divergence conditions up to discretization error. The strong
form dG discretization of equation (35) is: Find
$(\boldsymbol{H}_{h},\boldsymbol{E}_{h})\in P^{h}\times Q^{h}$ satisfying the
initial conditions (35e), such that

$\begin{split}\int_{\Omega}(\mu\boldsymbol{H}_{h,t}+\nabla\times\boldsymbol{E}_{h})\cdot\boldsymbol{G}_{h}\,d\boldsymbol{x}+\int_{\Omega}(\epsilon\boldsymbol{E}_{h,t}-\nabla\times\boldsymbol{H}_{h})\cdot\boldsymbol{F}_{h}\,d\boldsymbol{x}=\qquad\qquad\\\
\sum_{e}\\!\int_{\partial{\Omega^{e}}}\\!\\!\\!\\!\\!-\left(\boldsymbol{n}^{-}\\!\times\\!(\boldsymbol{E}_{h}^{\dagger}-\boldsymbol{E}_{h}^{-})\right)\cdot\boldsymbol{G}_{h}\,d\boldsymbol{x}+\sum_{e}\\!\int_{\partial{\Omega^{e}}}\\!\\!\\!\left(\boldsymbol{n}^{-}\times(\boldsymbol{H}_{h}^{\dagger}-\boldsymbol{H}_{h}^{-})\right)\cdot\boldsymbol{F}_{h}\,d\boldsymbol{x}\end{split}$
(36)

for all $(\boldsymbol{G}_{h},\boldsymbol{F}_{h})\in P^{h}\times Q^{h}$, and
for all $t\in(0,T)$. The upwind numerical flux states are given such that

$\displaystyle\boldsymbol{n}^{-}\times(\boldsymbol{E}_{h}^{\dagger}-\boldsymbol{E}_{h}^{-})$
$\displaystyle=-\frac{1}{2\\{\\!\\!\\{Y\\}\\!\\!\\}}\boldsymbol{n}^{-}\times(Y^{+}[\boldsymbol{E}_{h}]+\boldsymbol{n}^{-}\times[\boldsymbol{H}_{h}]),$
(37a)
$\displaystyle\boldsymbol{n}^{-}\times(\boldsymbol{H}_{h}^{\dagger}-\boldsymbol{H}_{h}^{-})$
$\displaystyle=-\frac{1}{2\\{\\!\\!\\{Z\\}\\!\\!\\}}\boldsymbol{n}^{-}\times(Z^{+}[\boldsymbol{H}_{h}]-\boldsymbol{n}^{-}\times[\boldsymbol{E}_{h}]).$
(37b)

The boundary conditions (35f) are imposed via the upwind numerical flux by
setting exterior values on the boundary of the domain $\Gamma$ such that

$\displaystyle\boldsymbol{H}_{h}^{+}$
$\displaystyle=-\boldsymbol{H}_{h}^{-}+2\boldsymbol{J}_{s},$ (38a)
$\displaystyle\boldsymbol{E}_{h}^{+}$ $\displaystyle=\boldsymbol{E}_{h}^{-},$
(38b)

with the continuously extended material parameters $Y^{+}=Y^{-}$ and
$Z^{+}=Z^{-}$. Using the upwind numerical flux implicitly means that the
boundary conditions are only set on the incoming characteristics.

Next, we compute the discrete adjoint equation and the gradient with respect
to the boundary current density $\boldsymbol{J}_{s}$ for an objective
functional of the form

$\tilde{\mathcal{J}}(\boldsymbol{J}_{s},\boldsymbol{E}_{h})=\int_{0}^{T}\\!\\!\int_{\Omega}j_{\Omega}(\boldsymbol{E}_{h})\,d\boldsymbol{x}\,dt+\int_{0}^{T}\\!\\!\int_{\Gamma}r(\boldsymbol{J}_{s})\,d\boldsymbol{x}\,dt,$
(39)

with differentiable functions $j_{\Omega}:\Omega\to\mathbb{R}$ and
$r:\Gamma\to\mathbb{R}$. We introduce the Lagrangian function, and derive the
adjoint equation by imposing that all variations of the Lagrangian with
respect to $\boldsymbol{H}_{h}$ and $\boldsymbol{E}_{h}$ must vanish. After
integration by parts in time $t$, this results in the following adjoint
equation: Find $(\boldsymbol{G}_{h},\boldsymbol{F}_{h})\in P^{h}\times Q^{h}$
such that

$\begin{split}\int_{\Omega}\mu\boldsymbol{G}_{h,t}\cdot\tilde{}\boldsymbol{H}_{h}+\boldsymbol{F}_{h}\cdot(\nabla\times\tilde{}\boldsymbol{H}_{h})\,d\boldsymbol{x}+\int_{\Omega}\epsilon\boldsymbol{F}_{h,t}\cdot\tilde{}\boldsymbol{E}_{h}-\boldsymbol{G}_{h}\cdot(\nabla\times\tilde{}\boldsymbol{E}_{h})\,d\boldsymbol{x}=\qquad\qquad\\\
\sum_{e}\\!\int_{\partial{\Omega^{e}}}\\!\\!\\!\\!-\left(\boldsymbol{n}^{-}\times\boldsymbol{F}_{h}^{\dagger}\right)\cdot\tilde{}\boldsymbol{H}_{h}\,d\boldsymbol{x}+\sum_{e}\\!\int_{\partial{\Omega^{e}}}\left(\boldsymbol{n}^{-}\times\boldsymbol{G}_{h}^{\dagger}\right)\cdot\tilde{}\boldsymbol{E}_{h}\,d\boldsymbol{x}-\int_{\Omega}j_{\Omega}^{\prime}(\boldsymbol{E}_{h})\cdot\tilde{}\boldsymbol{E}_{h}\,d\boldsymbol{x}\end{split}$
(40)

for all $(\tilde{}\boldsymbol{H}_{h},\tilde{}\boldsymbol{E}_{h})\in
P^{h}\times Q^{h}$ and all $t\in(0,T)$, with the final time conditions

$\boldsymbol{G}_{h}(\boldsymbol{x},T)=0,\>\boldsymbol{F}_{h}(\boldsymbol{x},T)=0\quad\text{on}\
\Omega,$ (41)

and the adjoint numerical flux states

$\displaystyle\boldsymbol{F}_{h}^{\dagger}$
$\displaystyle=\frac{\\{\\!\\!\\{Y\boldsymbol{F}_{h}\\}\\!\\!\\}}{\\{\\!\\!\\{Y\\}\\!\\!\\}}+\frac{1}{2\\{\\!\\!\\{Y\\}\\!\\!\\}}\left(\boldsymbol{n}^{-}\times[\boldsymbol{G}_{h}]\right),$
(42a) $\displaystyle\boldsymbol{G}_{h}^{\dagger}$
$\displaystyle=\frac{\\{\\!\\!\\{Z\boldsymbol{G}_{h}\\}\\!\\!\\}}{\\{\\!\\!\\{Z\\}\\!\\!\\}}-\frac{1}{2\\{\\!\\!\\{Z\\}\\!\\!\\}}\left(\boldsymbol{n}^{-}\times[\boldsymbol{F}_{h}]\right),$
(42b)

with exterior values on the boundary $\Gamma$ given by
$\boldsymbol{G}_{h}^{+}=-\boldsymbol{G}_{h}^{-}$ and
$\boldsymbol{F}_{h}^{+}=\boldsymbol{F}_{h}^{-}$. These exterior states enforce
the continuous adjoint boundary condition

$\boldsymbol{n}\times\boldsymbol{G}=0\quad\text{on}\ \Gamma.$

To compare with the numerical flux (37) of the state equation, we rewrite the
adjoint numerical flux states as

$\displaystyle\boldsymbol{n}^{-}\times(\boldsymbol{F}_{h}^{\dagger}-\boldsymbol{F}_{h}^{-})$
$\displaystyle=-\frac{1}{2\\{\\!\\!\\{Y\\}\\!\\!\\}}\boldsymbol{n}^{-}\times(Y^{+}[\boldsymbol{F}_{h}]-\boldsymbol{n}^{-}\times[\boldsymbol{G}_{h}]),$
$\displaystyle\boldsymbol{n}^{-}\times(\boldsymbol{G}_{h}^{\dagger}-\boldsymbol{G}_{h}^{-})$
$\displaystyle=-\frac{1}{2\\{\\!\\!\\{Z\\}\\!\\!\\}}\boldsymbol{n}^{-}\times(Z^{+}[\boldsymbol{G}_{h}]+\boldsymbol{n}^{-}\times[\boldsymbol{F}_{h}]).$

Note that, even with discontinuities in the material parameters, (40) is the
weak form of the dG discretization for a Maxwell’s system solved backwards in
time. As in the acoustic example, the adjoint numerical flux states (42) come
from the downwind flux.

Differentiating the Lagrangian with respect to the boundary current
$\boldsymbol{J}_{s}$ yields an equation for the derivative in direction
$\tilde{}\boldsymbol{J}_{s}$, namely

${\mathcal{J}}^{\prime}(\boldsymbol{J}_{s})(\tilde{}\boldsymbol{J}_{s})=\int_{0}^{T}\\!\\!\\!\int_{\Gamma}r^{\prime}(\boldsymbol{J}_{s})\cdot\tilde{}\boldsymbol{J}_{s}\,d\boldsymbol{x}+\int_{0}^{T}\\!\\!\int_{\Gamma}\boldsymbol{g}\cdot\tilde{}\boldsymbol{J}_{s}\,d\boldsymbol{x}\,dt$
(44a) with
$\boldsymbol{g}=\boldsymbol{n}^{-}\times\boldsymbol{F}_{h}+\frac{1}{Y}\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times\boldsymbol{G}_{h}\right)\right),$
(44b)

where $(\boldsymbol{G}_{h},\boldsymbol{F}_{h})$ is the solution of the
discrete adjoint equation (40).

###### Remark 5

Since the boundary force $\boldsymbol{J}_{s}$ enters linearly in Maxwell’s
equation, the gradient expression (44) does not involve contributions from the
element boundaries as for the advection and the acoustic wave example.

###  3.4 Coupled elastic-acoustic wave equation

Finally, we present expressions for the derivatives with respect to the
primary and secondary wave speeds in the coupled acoustic-elastic wave
equation. This section generalizes Section 3.2 to the coupled acoustic-elastic
wave equation. We derive derivative expressions with respect to both wave
speeds, and allow for discontinuous wave speeds across elements. We only
present a condensed form of the derivations, and verify our results for the
gradient numerically in Section 4.

The coupled linear elastic-acoustic wave equation for isotropic material
written in first-order velocity strain form is given as

$\displaystyle\boldsymbol{E}_{t}$
$\displaystyle=\frac{1}{2}\left(\nabla\boldsymbol{v}+\nabla\boldsymbol{v}^{T}\right)$
$\displaystyle\qquad\text{\ on\ }\Omega\times(0,T)$ (45a)
$\displaystyle\rho\boldsymbol{v}_{t}$
$\displaystyle=\nabla\cdot\left(\lambda\operatorname{tr}(\boldsymbol{E})\boldsymbol{I}+2\mu\boldsymbol{E}\right)+\rho\boldsymbol{f}$
$\displaystyle\qquad\text{\ on\ }\Omega\times(0,T)$ (45b) where
$\boldsymbol{E}$ is the strain tensor, $\boldsymbol{v}$ is the displacement
velocity, $\boldsymbol{I}$ is the identity tensor, $\rho=\rho(\boldsymbol{x})$
is the mass density, $\boldsymbol{f}$ is a body force per unit mass, and
$\lambda=\lambda(\boldsymbol{x})$ and $\mu=\mu(\boldsymbol{x})$ are the Lamé
parameters. In addition to the conditions (45a) and (45b) on the body
$\Omega$, we assume the initial conditions
$\displaystyle\boldsymbol{v}(0,\boldsymbol{x})=\boldsymbol{v}_{0}(\boldsymbol{x}),\>\boldsymbol{E}(0,\boldsymbol{x})=\boldsymbol{E}_{0}(\boldsymbol{x}),\quad\text{\
for\ }\boldsymbol{x}\in\Omega,$ (45c) and the boundary conditions
$\displaystyle\SS(\boldsymbol{x},t)\boldsymbol{n}=\boldsymbol{t}^{\text{bc}}(t)\qquad\text{on}\
\Gamma.$ (45d) Here, $\boldsymbol{t}^{\text{bc}}$ is the traction on the
boundary of the body.

The stress tensor $\SS$ is related to the strain through the constitutive
relation (here, $\boldsymbol{\mathsf{C}}$ is the forth-order constitutive
tensor):

$\SS=\boldsymbol{\mathsf{C}}\boldsymbol{E}=\lambda\operatorname{tr}(\boldsymbol{E})\boldsymbol{I}+2\mu\boldsymbol{E},$
(46)

where $\operatorname{tr}(\cdot)$ is the trace operator. There are also
boundary conditions at material interfaces. For an elastic-elastic interface
$\Gamma^{\text{ee}}$ the boundary conditions are

$\displaystyle\boldsymbol{v}^{+}=\boldsymbol{v}^{-},\quad\SS^{+}\boldsymbol{n}^{-}$
$\displaystyle=\SS^{-}\boldsymbol{n}^{-}\qquad\text{on}\ \Gamma^{\text{ee}}$

and for acoustic-elastic and acoustic-acoustic interfaces
$\Gamma^{\text{ae}}$, the boundary conditions are

$\displaystyle\boldsymbol{n}\cdot\boldsymbol{v}^{+}=\boldsymbol{n}\cdot\boldsymbol{v}^{-},\quad\SS^{+}\boldsymbol{n}^{-}=\SS^{-}\boldsymbol{n}^{-}\qquad\text{on}\
\Gamma^{\text{ae}}.$

The strong form dG discretization of equation (45) is: Find
$(\boldsymbol{E}_{h},\boldsymbol{v}_{h})\in P^{h}\times Q^{h}$ such that

$\int_{\Omega}\boldsymbol{E}_{h,t}:\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}\,d\boldsymbol{x}+\int_{\Omega}\rho\boldsymbol{v}_{h,t}\cdot\boldsymbol{w}_{h}\,d\boldsymbol{x}-\int_{\Omega}\operatorname{sym}(\nabla\boldsymbol{v}_{h}):\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}\,d\boldsymbol{x}\\\
-\int_{\Omega}\left(\nabla\cdot(\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h})+\boldsymbol{f}\right)\cdot\boldsymbol{w}_{h}\,d\boldsymbol{x}=\sum_{e}\int_{\partial{\Omega^{e}}}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{v}_{h}^{\dagger}-\boldsymbol{v}_{h}^{-}\right)\right):\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}\,d\boldsymbol{x}\\\
+\sum_{e}\int_{\partial{\Omega^{e}}}\left(\left({(\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h})}^{\dagger}-{(\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h})}^{-}\right)\boldsymbol{n}^{-}\right)\cdot\boldsymbol{w}_{h}\,d\boldsymbol{x}$
(49)

for all $(\boldsymbol{H}_{h},\boldsymbol{w}_{h})\in P^{h}\times Q^{h}$ where
$\operatorname{sym}$ is the mapping to get the symmetric part of a tensor,
i.e.,
$\operatorname{sym}(\boldsymbol{A})=\frac{1}{2}\left(\boldsymbol{A}+\boldsymbol{A}^{T}\right)$.
Note that the constitutive tensor $\boldsymbol{\mathsf{C}}$ is used in the
inner product for the weak form. Here, the upwind states are given such that

$\displaystyle\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{v}_{h}^{\dagger}-\boldsymbol{v}_{h}^{-}\right)\right)$
$\displaystyle={-k_{0}\left(\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!]+\rho^{+}c_{p}^{+}[\\![\boldsymbol{v}_{h}]\\!]\right)}\left(\boldsymbol{n}^{-}\otimes\boldsymbol{n}^{-}\right)$
$\displaystyle\quad+k_{1}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!]\right)\right)\right)$
$\displaystyle\quad+k_{1}\rho^{+}c_{s}^{+}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\boldsymbol{v}_{h}]\right)\right)\right),$
(50a)
$\displaystyle\left({(\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h})}^{\dagger}-{(\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h})}^{-}\right)\boldsymbol{n}^{-}$
$\displaystyle={-k_{0}\left(\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!]+\rho^{+}c_{p}^{+}[\\![\boldsymbol{v}_{h}]\\!]\right)}\rho^{-}c_{p}^{-}\boldsymbol{n}^{-}$
$\displaystyle\quad+k_{1}\rho^{-}c_{s}^{-}\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!]\right)$
$\displaystyle\quad+k_{1}\rho^{+}c_{s}^{+}\rho^{-}c_{s}^{-}\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\boldsymbol{v}_{h}]\right),$
(50b) with $k_{0}=1/(\rho^{-}c_{p}^{-}+\rho^{+}c_{p}^{+})$ and $\displaystyle
k_{1}$
$\displaystyle=\begin{dcases}\frac{1}{\rho^{-}c_{s}^{-}+\rho^{+}c_{s}^{+}}&\text{when}\
\mu^{-}\neq 0,\\\ 0&\text{when}\ \mu^{-}=0,\end{dcases}$

where $c_{p}:=\sqrt{(\lambda+2\mu)/\rho}$ is the primary wave speed and
$c_{s}:=\sqrt{\mu/\rho}$ is the secondary wave speed. The traction boundary
conditions are imposed through the upwind numerical flux by setting exterior
values on $\Gamma$ to

$\displaystyle\boldsymbol{v}_{h}^{+}$ $\displaystyle=\boldsymbol{v}_{h}^{-},$
$\displaystyle\boldsymbol{\mathsf{C}}^{+}\boldsymbol{E}_{h}^{+}\boldsymbol{n}^{+}$
$\displaystyle=\begin{cases}-2\boldsymbol{t}^{\text{bc}}+\boldsymbol{\mathsf{C}}^{-}\boldsymbol{E}_{h}^{-}\boldsymbol{n}^{-}&\
\text{if $\mu^{-}\not=0$},\\\
-2\left(\boldsymbol{n}^{-}\cdot\left(\boldsymbol{t}^{\text{bc}}-\boldsymbol{\mathsf{C}}^{-}\boldsymbol{E}_{h}^{-}\boldsymbol{n}^{-}\right)\right)\boldsymbol{n}^{-}-\boldsymbol{\mathsf{C}}^{-}\boldsymbol{E}_{h}^{-}\boldsymbol{n}^{-}&\
\text{if $\mu^{+}=0$},\end{cases}$

with the continuously extended material parameters $\rho^{+}=\rho^{-}$,
$\mu^{+}=\mu^{-}$, and $\lambda^{+}=\lambda^{-}$.

We assume a cost function that depends on the primary and secondary wave
speeds $c_{p}$ and $c_{s}$ through the solution $\boldsymbol{v}_{h}$ of (49)

${\mathcal{J}}(c_{p},c_{s})=\int_{0}^{T}\\!\\!\\!\int_{\Omega}j_{\Omega}(\boldsymbol{v}_{h})\,d\boldsymbol{x}\,dt+\int_{\Omega}r_{p}(c_{p})\,d\boldsymbol{x}+\int_{\Omega}r_{s}(c_{s})\,d\boldsymbol{x}.$
(51)

By using a sum of spatial Dirac delta distributions in $j_{\Omega}(\cdot)$,
this can include seismogram data, as common in seismic inversion. Using the
Lagrangian function and integration by parts in time, we obtain the following
adjoint equation: Find $(\boldsymbol{H}_{h},\boldsymbol{w}_{h})\in P^{h}\times
Q^{h}$ such that

$\int_{\Omega}-\boldsymbol{H}_{h,t}:\boldsymbol{\mathsf{C}}\tilde{\boldsymbol{E}}_{h}\,d\boldsymbol{x}-\int_{\Omega}\rho\boldsymbol{w}_{h,t}\cdot\tilde{\boldsymbol{v}}_{h}\,d\boldsymbol{x}-\int_{\Omega}\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}:\operatorname{sym}(\nabla\tilde{\boldsymbol{v}}_{h})\,d\boldsymbol{x}\\\
-\int_{\Omega}\boldsymbol{w}_{h}\cdot\left(\nabla\cdot(\boldsymbol{\mathsf{C}}\tilde{\boldsymbol{E}}_{h})\right)\,d\boldsymbol{x}=-\sum_{e}\int_{\partial{\Omega^{e}}}\operatorname{sym}(\boldsymbol{n}^{-}\otimes\boldsymbol{w}_{h}^{\dagger}):\boldsymbol{\mathsf{C}}\tilde{\boldsymbol{E}}_{h}\,d\boldsymbol{x}\\\
-\sum_{e}\int_{\partial{\Omega^{e}}}\left({(\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h})}^{\dagger}\boldsymbol{n}^{-}\right)\cdot\tilde{\boldsymbol{v}}_{h}\,d\boldsymbol{x},-\int_{\Omega}j_{\Omega}(\boldsymbol{v}_{h})\cdot\tilde{\boldsymbol{v}}_{h}$
(52)

for all $(\tilde{\boldsymbol{E}}_{h},\tilde{\boldsymbol{v}}_{h})\in
P^{h}\times Q^{h}$ in with final conditions $\boldsymbol{w}_{h}(T)=0$ and
$\boldsymbol{H}_{h}(T)=0$ the fluxes are given by

$\displaystyle\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\boldsymbol{w}_{h}^{\dagger}\right)$
$\displaystyle=k_{0}\left(\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]+\boldsymbol{n}^{-}\cdot\left(2\\{\\!\\!\\{\rho
c_{p}\boldsymbol{w}_{h}\\}\\!\\!\\}\right)\right)\left(\boldsymbol{n}^{-}\otimes\boldsymbol{n}^{-}\right)$
$\displaystyle\quad-
k_{1}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]\right)\right)\right)$
$\displaystyle\quad-
k_{1}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times\left(2\\{\\!\\!\\{\rho
c_{s}\boldsymbol{w}_{h}\\}\\!\\!\\}\right)\right)\right)\right),$
$\displaystyle{\left(\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}\right)}^{\dagger}\boldsymbol{n}^{-}$
$\displaystyle=k_{0}\left(\boldsymbol{n}^{-}\cdot\left(\left(\rho^{+}c_{p}^{+}\boldsymbol{\mathsf{C}}^{-}\boldsymbol{H}_{h}^{-}+\rho^{-}c_{p}^{-}\boldsymbol{\mathsf{C}}^{+}\boldsymbol{H}_{h}^{+}\right)\boldsymbol{n}^{-}\right)+\rho^{-}c_{p}^{-}\rho^{+}c_{p}^{+}[\\![\boldsymbol{w}_{h}]\\!]\right)\boldsymbol{n}^{-}$
$\displaystyle\quad-
k_{1}\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times\left(\left(\rho^{+}c_{s}^{+}\boldsymbol{\mathsf{C}}^{-}\boldsymbol{H}_{h}^{-}+\rho^{-}c_{s}^{-}\boldsymbol{\mathsf{C}}^{+}\boldsymbol{H}_{h}^{+}\right)\boldsymbol{n}^{-}\right)\right)$
$\displaystyle\quad-
k_{1}\rho^{-}c_{s}^{-}\rho^{+}c_{s}^{+}\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\boldsymbol{w}_{h}]\right).$

We can rewrite this into a form similar to the upwind states of the state
equation (50) as

$\displaystyle\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{w}_{h}^{\dagger}-\boldsymbol{w}_{h}^{-}\right)\right)$
$\displaystyle=k_{0}\left(\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]-\rho^{+}c_{p}^{+}[\\![\boldsymbol{w}_{h}]\\!]\right)\left(\boldsymbol{n}^{-}\otimes\boldsymbol{n}^{-}\right)$
$\displaystyle\quad-
k_{1}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]\right)\right)\right)$
$\displaystyle\quad+k_{1}\rho^{+}c_{s}^{+}\operatorname{sym}\left(\boldsymbol{n}^{-}\otimes\left(\boldsymbol{n}^{-}\times\left(\boldsymbol{n}^{-}\times[\boldsymbol{w}_{h}]\right)\right)\right),$
$\displaystyle\left({(\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h})}^{\dagger}-{(\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h})}^{-}\right)\boldsymbol{n}^{-}$
$\displaystyle=k_{0}\left(-\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]+\rho^{+}c_{p}^{+}[\\![\boldsymbol{w}_{h}]\\!]\right)\rho^{-}c_{p}^{-}\boldsymbol{n}^{-}$
$\displaystyle\quad+k_{1}\rho^{-}c_{s}^{-}\boldsymbol{n}^{-}\times(\boldsymbol{n}^{-}\times[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!])$
$\displaystyle\quad-
k_{1}\rho^{-}c_{s}^{-}\rho^{+}c_{s}^{+}\boldsymbol{n}^{-}\times(\boldsymbol{n}^{-}\times[\boldsymbol{w}_{h}]).$

Here, the adjoint boundary conditions are imposed through the adjoint
numerical flux by setting exterior values on $\Gamma$ to

$\displaystyle\boldsymbol{w}_{h}^{+}$ $\displaystyle=\boldsymbol{w}_{h}^{-},$
$\displaystyle\boldsymbol{\mathsf{C}}^{+}\boldsymbol{H}_{h}^{+}\boldsymbol{n}^{+}$
$\displaystyle=\begin{cases}\phantom{-}\boldsymbol{\mathsf{C}}^{-}\boldsymbol{H}_{h}^{-}\boldsymbol{n}^{-}&\
\text{if $\mu^{-}\not=0$},\\\
-\boldsymbol{\mathsf{C}}^{-}\boldsymbol{H}_{h}^{-}\boldsymbol{n}^{-}+2\left(\boldsymbol{n}^{-}\cdot\left(\boldsymbol{\mathsf{C}}^{-}\boldsymbol{H}_{h}^{-}\boldsymbol{n}^{-}\right)\right)\boldsymbol{n}^{-}&\
\text{if $\mu^{-}=0$},\end{cases}$

with the continuously extended material parameters $\rho^{+}=\rho^{-}$,
$\mu^{+}=\mu^{-}$, and $\lambda^{+}=\lambda^{-}$.

We assume a discretization of $c_{p}$ and $c_{s}$ and a numerical quadrature
rule such that the state equation can be used to simplify the expression for
the gradient; see the discussion in Section 3.2. The discrete gradient with
respect to $c_{p}$ is then

${\mathcal{J}}_{c_{p}}(c_{p},c_{s})(\tilde{c}_{p})=\int_{\Omega}r_{p}^{\prime}(c_{p})\tilde{c}_{p}\,d\boldsymbol{x}+\int_{0}^{T}\\!\\!\int_{\Omega}G_{p}\tilde{c}_{p}\,d\boldsymbol{x}\,dt+\sum_{\partial{\Omega^{e}}}\int_{0}^{T}\\!\\!\int_{\partial{\Omega^{e}}}g_{p}\tilde{c}^{-}_{p}\,d\boldsymbol{x}\,dt,$
(53a) where $G_{p}$ is defined on each element ${\Omega^{e}}$, and $g_{p}$ for
each element boundary as follows:
$\begin{split}G_{p}\tilde{c}_{p}&=-2\left(\nabla\left(\rho{c_{p}}{\tilde{c}_{p}}\operatorname{tr}(\boldsymbol{E}_{h})\right)\right)\cdot\boldsymbol{w}_{h},\\\
g_{p}&=-k_{0}^{2}\rho^{-}\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!]\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]+k_{0}^{2}\rho^{-}{\left(\rho^{+}{c_{p}}^{+}\right)}^{2}[\\![\boldsymbol{v}_{h}]\\!][\\![\boldsymbol{w}_{h}]\\!]\\\
&\quad+k_{0}^{2}\rho^{-}\rho^{+}{c_{p}}^{+}\left(\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!][\\![\boldsymbol{w}_{h}]\\!]-[\\![\boldsymbol{v}_{h}]\\!]\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]\right)\\\
&\quad+2k_{0}\rho^{-}{c_{p}}^{-}\operatorname{tr}(\boldsymbol{E}_{h}^{-})\left(\boldsymbol{n}^{-}\cdot[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]-\rho^{+}{c_{p}}^{+}[\\![\boldsymbol{w}_{h}]\\!]\right)\\\
&\quad+2\rho^{-}c_{p}^{-}\operatorname{tr}(\boldsymbol{E}_{h}^{-})\boldsymbol{w}_{h}^{-}\cdot\boldsymbol{n}^{-},\end{split}$
(53b)

where $(\boldsymbol{v}_{h},\boldsymbol{E}_{h})$ is the solution of the state
equation (49) and $(\boldsymbol{w}_{h},\boldsymbol{H}_{h})$ is the solution of
the adjoint equation (52). The discrete gradient of $\mathcal{J}$ with respect
to $c_{s}$ is

${\mathcal{J}}_{c_{s}}(c_{p},c_{s})(\tilde{c}_{s})=\int_{\Omega}r_{s}^{\prime}(c_{s})\tilde{c}_{s}\,d\boldsymbol{x}+\int_{0}^{T}\\!\\!\int_{\Omega}G_{s}\tilde{c}_{s}\,d\boldsymbol{x}\,dt+\sum_{\partial{\Omega^{e}}}\int_{0}^{T}\\!\\!\int_{\partial{\Omega^{e}}}g_{s}\tilde{c}^{-}_{s}\,d\boldsymbol{x}\,dt,$
(54a) where $G_{s}$ is defined on each element ${\Omega^{e}}$, and $g_{s}$ for
each element boundary as follows:
$\begin{split}G_{s}\tilde{c}_{s}&=-4\left(\nabla\cdot\left(\rho
c_{s}\tilde{c}_{s}\left(\boldsymbol{E}_{h}-\operatorname{tr}(\boldsymbol{E}_{h})\boldsymbol{I}\right)\right)\right)\cdot\boldsymbol{w}_{h},\\\
g_{s}&=-k_{1}^{2}\rho^{-}\\!\left(\boldsymbol{n}^{-}\\!\times\\!\left(\boldsymbol{n}^{-}\\!\times{[\\![\boldsymbol{\mathsf{C}}\boldsymbol{E}_{h}]\\!]}\right)\right)\cdot\left(\boldsymbol{n}^{-}\\!\times\\!\left(\boldsymbol{n}^{-}\\!\times{\big{(}[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]-\rho^{+}c_{s}^{+}[\boldsymbol{w}_{h}]\big{)}}\right)\right)\\\
&\quad-
k_{1}^{2}\rho^{-}\rho^{+}{c_{s}}^{+}\\!\left(\boldsymbol{n}^{-}\\!\times\\!\left(\boldsymbol{n}^{-}\\!\times{[\boldsymbol{v}_{h}]}\right)\right)\cdot\left(\boldsymbol{n}^{-}\\!\times\\!\left(\boldsymbol{n}^{-}\\!\times{\big{(}[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]-\rho^{+}c_{s}^{+}[\boldsymbol{w}_{h}]\big{)}}\right)\right)\\\
&\quad+4k_{0}\rho^{-}{c_{s}}^{-}\\!\left(\boldsymbol{n}^{-}\cdot\boldsymbol{E}_{h}^{-}\boldsymbol{n}^{-}-\operatorname{tr}(\boldsymbol{E}_{h}^{-})\right)\left(\boldsymbol{n}^{-}\cdot\left([\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]-\rho^{+}{c_{p}}^{+}[\boldsymbol{w}_{h}]\right)\right)\\\
&\quad+4k_{1}\rho^{-}{c_{s}}^{-}\\!\left(\boldsymbol{n}^{-}\\!\times\\!\left(\boldsymbol{n}^{-}\\!\times{\big{(}\boldsymbol{E}_{h}^{-}\boldsymbol{n}^{-}\big{)}}\right)\right)\\!\cdot\\!\left(\boldsymbol{n}^{-}\\!\times\\!\left(\boldsymbol{n}^{-}\\!\times{\big{(}[\\![\boldsymbol{\mathsf{C}}\boldsymbol{H}_{h}]\\!]\\!-\\!\rho^{+}{c_{s}}^{+}[\boldsymbol{w}_{h}]\big{)}}\right)\right)\\\
&\quad+4\rho^{-}c_{s}^{-}\\!\left(\left(\boldsymbol{E}_{h}^{-}-\operatorname{tr}(\boldsymbol{E}_{h}^{-})\boldsymbol{I}\right)\boldsymbol{n}^{-}\right)\cdot\boldsymbol{w}_{h}^{-},\end{split}$
(54b)

where $(\boldsymbol{v}_{h},\boldsymbol{E}_{h})$ is the solution of the state
equation (49) and $(\boldsymbol{w}_{h},\boldsymbol{H}_{h})$ is the solution of
the adjoint equation (52). Note that since we allow for discontinuous wave
speeds $c_{p}$ and $c_{s}$ (and perturbations $\tilde{c}_{p}$ and
$\tilde{c}_{s}$), the boundary contributions to the gradients, i.e., the last
terms in (53a) and (54a) are written as sums of integrals over individual
element boundaries, i.e., each boundary face appears twice in the overall sum.
This differs from the previous examples, where we assumed continuous
parameters and thus combined contributions from adjacent elements to shared
faces $\mathsf{f}$.

###### Remark 6

The expressions for the $c_{p}$-gradient (53) reduce to the result found for
the acoustic equation (32) and (34) for continuous parameter fields
$\rho,c_{p},c_{s}$ and continuous parameter perturbations $\tilde{c}_{p}$. To
verify this, one adds contributions from adjacent elements to common
boundaries, and the terms in (53b) combine or cancel.

###### Remark 7

Above, we have derived expressions for the derivatives with respect to the
primary and secondary wave speeds. If, instead of $c_{p}$ and $c_{s}$,
derivatives with respect to an alternative pair of parameters in the stress
tensor—such as the Lamé parameters, or Poisson’s ratio and Young’s modulus—are
derived, the adjoint equations remain unchanged, but the expressions for the
derivatives change according to the chain rule.

##  4 Numerical verification for coupled elastic-acoustic wave propagation

Here, we numerically verify the expressions for the discrete gradients with
respect to the wave speeds for the elastic-acoustic wave problem derived in
Section 3.4. For this purpose, we compare directional finite differences with
directional gradients based on the discrete gradient. To emphasize the
correctness of the discrete gradient, we use coarse meshes in these
comparisons, which underresolve the wave fields. As test problem, we use the
Snell law example from Section 6.2 in WilcoxStadlerBursteddeEtAl10  with the
material parameters and the wave incident angle specified there. For our
tests, we use the simple distributed objective function

$\mathcal{J}(c_{p},c_{s}):=\int_{0}^{T}\\!\\!\int_{\Omega}\boldsymbol{v}_{h}\cdot\boldsymbol{v}_{h}\,d\boldsymbol{x}\,dt$

The discretization of the wave equation follows WilcoxStadlerBursteddeEtAl10 ,
i.e., we use spectral elements based on Gauss-Lobatto-Lagrange (GLL) points on
hexahedral meshes. The use of GLL quadrature results in underintegration even
if the elements are images of the reference element under an affine
transformation. In Figure 2, we summarize results for the directional
derivatives in the direction $\tilde{c}_{p}:=\sin(\pi x)\cos(\pi y)\cos(\pi
z)$; we compare the finite difference directional derivatives

$d_{\epsilon}^{\text{fd}}:=\frac{\mathcal{J}(c_{p}+\epsilon\tilde{c}_{p})-\mathcal{J}(c_{p})}{\epsilon}$
(55)

with the directional derivatives $d^{\text{di}}$ and $d^{\text{co}}$ defined
by

$d^{\text{di}}:=\mathcal{J}_{c_{p}}(c_{p},c_{s})(\tilde{c}_{p}),\quad
d^{\text{co}}:=\mathcal{J}^{\text{cont}}_{c_{p}}(c_{p},c_{s})(\tilde{c}_{p}),$
(56)

where $\mathcal{J}_{c_{p}}(c_{p},c_{s})$ denotes the discrete gradient (53),
and $\mathcal{J}^{\text{co}}_{c_{p}}(c_{p},c_{s})$ denotes the gradient
obtained when neglecting the jump term in the boundary contributions $g_{p}$
in (53b). These jump terms are likely to be neglected if the continuous
gradient expressions are discretized instead of following a fully discrete
approach. The resulting error is of the order of the discretization error and
thus vanishes as the discrete solutions converge. However, this error can be
significant on coarse meshes, on which the wave solution is not well resolved.

$1$$2$$3$$0.9$$1$$1.1$$1.2$$1.3$mesh leveldirectional
derivative$d^{\text{fd}}_{\epsilon},\epsilon=10^{-3}$$d^{\text{fd}}_{\epsilon},\epsilon=10^{-4}$$d^{\text{fd}}_{\epsilon},\epsilon=10^{-5}$disc.
grad. $d^{\text{di}}$cont. grad. $d^{\text{co}}$

mesh level 1
---
$d^{\text{fd}}_{\epsilon},\epsilon=10^{-2}$ |  1.489022
$d^{\text{fd}}_{\epsilon},\epsilon=10^{-3}$ |  1.244358
$d^{\text{fd}}_{\epsilon},\epsilon=10^{-4}$ |  1.219841
$d^{\text{fd}}_{\epsilon},\epsilon=10^{-5}$ |  1.217389
$d^{\text{fd}}_{\epsilon},\epsilon=10^{-6}$ |  1.217143
$d^{\text{fd}}_{\epsilon},\epsilon=10^{-7}$ |  1.217118
$d^{\text{di}}$ | 1.217117

Figure 2: Directional derivatives computed using one-sided finite differences
(55), and the discrete and the continuous gradients (56). Left: Results on
mesh levels 1,2 and 3 corresponding to meshes with 16, 128 and 1024 finite
elements with polynomial order $N=4$. The finite difference directional
derivatives $d^{\text{fd}}_{\epsilon}$ converge to the discrete gradient
$d^{\text{di}}$ as $\epsilon$ is reduced. Note that as the mesh level is
increased, the continuous gradient $d^{\text{co}}$ converges to
$d^{\text{di}}$. Right: Convergence of finite difference directional
derivative on the coarsest mesh. Digits for which the finite difference
gradient coincides with the discrete gradient are shown in bold.

Next, we study the accuracy of the discrete gradient for pointwise
perturbations to the wave speed. Since the same discontinuous basis functions
as for the wave equation are also used for the local wave speeds, a point
perturbation in $c_{p}^{-}$ or $c_{s}^{-}$ at an element boundary face results
in a globally discontinuous perturbation direction $\tilde{c}_{p}$ and
$\tilde{c}_{s}$. In Table 1, we present the discrete directional gradient
$d^{\text{di}}$ with finite difference directional gradients
$d^{\text{fd}}_{\epsilon}$ for unit vector perturbations of both wave speeds.
Compared to in the table in Figure 2, where the directional derivatives for
smooth perturbations are reported, pointwise perturbations of the wave speeds
$c_{p}$ or $c_{s}$ result in smaller changes in the cost functional, and
numerical roundoff influences the accuracy of finite difference directional
derivatives. As a consequence, fewer digits coincide between the finite
difference directional derivatives and the discrete gradients.

Table 1: Comparison of pointwise material gradients for Snell problem from (WilcoxStadlerBursteddeEtAl10, , Section 6.2). The derivatives $d^{\text{fd}}_{\epsilon}$ and $d^{\text{di}}$ with respect to the local wave speed (either $c_{p}$ or $c_{s}$) for points with coordinates $(x,y,z)$ are reported. We use the final time $T=1$ and spectral elements of polynomial order $N=6$ in space. The meshes for level 1 and 2 consist of 16 and 128 finite elements, respectively. Digits where the finite difference approximation coincides with the discrete gradient are shown in bold. mesh level | $(x,y,z)$ | pert. | $d^{\text{di}}$ | $d^{\text{fd}}_{\epsilon}$
---|---|---|---|---
$\\#$tsteps |  | field |  | $\epsilon=10^{-3}$ | $\epsilon=10^{-4}$ | $\epsilon=10^{-5}$
1/101 | $(0,0,0)$ | $c_{p}$ | 1.8590e-4 |  1.8549e-4 |  1.8581e-4 |  1.8560e-4
2/202 | $(0,0,0)$ | $c_{p}$ | 2.2102e-5 |  2.2094e-5 |  2.2007e-5 |  2.1504e-5
1/101 | $(0,0,1)$ | $c_{s}$ | 1.1472e-5 |  1.1453e-5 |  1.1372e-5 |  1.0942e-5
1/101 | $(-0.5,-0.5,0.5)$ | $c_{s}$ | 2.8886e-3 |  2.8802e-3 |  2.8877e-3 |  2.8870e-3

##  5 Conclusions

Our study yields that the discretely exact adjoint PDE of a dG-discretized
linear hyperbolic equation is a proper dG discretization of the continuous
adjoint equation, provided an upwind flux is used. Thus, the adjoint PDE
converges at the same rate as the state equation. When integration by parts is
avoided to eliminate quadrature errors, a weak dG discretization of the state
PDE leads to a strong dG discretization of the adjoint PDE, and vice versa.
The expressions for the discretely exact gradient can contain contributions at
element faces and, hence, differ from a straightforward discretization of the
continuous gradient expression. These element face contributions are at the
order of the discretization order and are thus more significant for poorly
resolved state PDEs. We believe that these observations are relevant for
inverse problems and optimal control problems governed by hyperbolic PDEs
discretized by the discontinuous Galerkin method.

## Acknowledgments

We would like to thank Jeremy Kozdon and Gregor Gassner for fruitful
discussions and helpful comments, and Carsten Burstedde for his help with the
implementation of the numerical example presented in Section 4. Support for
this work was provided through the U.S. National Science Foundation (NSF)
grant CMMI-1028889, the Air Force Office of Scientific Research’s
Computational Mathematics program under the grant FA9550-12-1-0484, and
through the Mathematical Multifaceted Integrated Capability Centers (MMICCs)
effort within the Applied Mathematics activity of the U.S. Department of
Energy’s Advanced Scientific Computing Research program, under Award Number
DE-SC0009286. The views expressed in this document are those of the authors
and do not reflect the official policy or position of the Department of
Defense or the U.S. Government.

## References

  * (1) M. Alexe and A. Sandu, Space-time adaptive solution of inverse problems with the discrete adjoint method, Tech. Rep. 10-14, Virginia Tech, 2010\.
  * (2) R. Becker, D. Meidner, and B. Vexler, Efficient numerical solution of parabolic optimization problems by finite element methods, Optimization Methods Software, 22 (2007), pp. 813–833.
  * (3) A. Borzì and V. Schulz, Computational Optimization of Systems Governed by Partial Differential Equations, SIAM, 2012.
  * (4) M. Braack, Optimal control in fluid mechanics by finite elements with symmetric stabilization, SIAM Journal on Control and Optimization, 48 (2009), pp. 672–687.
  * (5) T. Bui-Thanh, C. Burstedde, O. Ghattas, J. Martin, G. Stadler, and L. C. Wilcox, Extreme-scale UQ for Bayesian inverse problems governed by PDEs, in SC12: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 2012.
  * (6) T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler, A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion, SIAM Journal on Scientific Computing, 35 (2013), pp. A2494–A2523.
  * (7) S. S. Collis and M. Heinkenschloss, Analysis of the streamline upwind/Petrov Galerkin method applied to the solution of optimal control problems, Tech. Rep. TR02–01, Department of Computational and Applied Mathematics, Rice University, Houston, TX 77005–1892, 2002.  http://www.caam.rice.edu/~heinken.
  * (8) S. S. Collis, C. C. Ober, and B. G. van Bloemen Waanders, Unstructured discontinuous Galerkin for seismic inversion, in Conference Paper, SEG International Exposition and 80th Annual Meeting, 2010.
  * (9) I. Epanomeritakis, V. Akçelik, O. Ghattas, and J. Bielak, A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion, Inverse Problems, 24 (2008), p. 034015 (26pp).
  * (10) K.-A. Feng, C.-H. Teng, and M.-H. Chen, A pseudospectral penalty scheme for 2D isotropic elastic wave computations, Journal of Scientific Computing, 33 (2007), pp. 313–348.
  * (11) A. Fichtner, Full seismic waveform modelling and inversion, Springer, 2011.
  * (12) M. Giles and N. Pierce, Adjoint equations in CFD: Duality, boundary conditions and solution behaviour, American Institute of Aeronautics and Astronautics, (1997).
  * (13) M. Giles and S. Ulbrich, Convergence of linearized and adjoint approximations for discontinuous solutions of conservation laws. Part 1: Linearized approximations and linearized output functionals, SIAM Journal on Numerical Analysis, 48 (2010), pp. 882–904.
  * (14) , Convergence of linearized and adjoint approximations for discontinuous solutions of conservation laws. Part 2: Adjoint approximations and extensions, SIAM Journal on Numerical Analysis, 48 (2010), pp. 905–921.
  * (15) A. Griewank and A. Walther, Algorithm 799: revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation, ACM Trans. Math. Softw., 26 (2000), pp. 19–45.
  * (16) A. Griewank and A. Walther, Evaluating derivatives: Principles and techniques of algorithmic differentiation, Society for Industrial Mathematics, second ed., 2008.
  * (17) M. D. Gunzburger, Perspectives in Flow Control and Optimization, SIAM, Philadelphia, 2003.
  * (18) B. Gustafsson, H.-O. Kreiss, and J. Oliger, Time Dependent Problems and Difference Methods, John Wiley & Sons, 1995.
  * (19) W. W. Hager, Runge-Kutta methods in optimal control and the transformed adjoint system, Numerische Mathematik, 87 (2000), pp. 247–282.
  * (20) K. Harriman, D. Gavaghan, and E. Süli, The importance of adjoint consistency in the approximation of linear functionals using the discontinuous Galerkin finite element method, Tech. Rep. NA04/18, Oxford University Computing Laboratory, 2004.
  * (21) R. Hartmann, Adjoint consistency analysis of discontinuous Galerkin discretizations, SIAM Journal on Numerical Analysis, 45 (2007), pp. 2671–2696.
  * (22) J. S. Hesthaven and T. Warburton, Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell’s equations, Journal of Computational Physics, 181 (2002), pp. 186–221.
  * (23) J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, vol. 54 of Texts in Applied Mathematics, Springer, 2008.
  * (24) M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich, Optimization with PDE Constraints, Springer, 2009.
  * (25) D. A. Kopriva, Implementing Spectral Methods for Partial Differential Equations, Springer, 2009.
  * (26) D. A. Kopriva and G. Gassner, On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods, Journal of Scientific Computing, 44 (2010), pp. 136–155.
  * (27) J. E. Lagnese, Exact boundary controllability of Maxwell’s equations in a general region, SIAM Journal on Control and Optimization, 27 (1989), pp. 374–388.
  * (28) V. Lekić and B. Romanowicz, Inferring upper-mantle structure by full waveform tomography with the spectral element method, Geophysical Journal International, 185 (2011), pp. 799–831.
  * (29) D. Leykekhman, Investigation of commutative properties of discontinuous Galerkin methods in PDE constrained optimal control problems, Journal of Scientific Computing, 53 (2012), pp. 483–511.
  * (30) J.-L. Lions, Control of Distributed Singular Systems, Dunod, Gauthier–Villars, Paris, 1985.
  * (31) S. Nicaise, Exact boundary controllability of Maxwell’s equations in heterogeneous media and an application to an inverse source problem, SIAM Journal on Control and Optimization, 38 (2000), pp. 1145–1170.
  * (32) T. A. Oliver and D. L. Darmofal, Analysis of dual consistency for discontinuous Galerkin discretizations of source terms, SIAM Journal on Numerical Analysis, 47 (2009), pp. 3507–3525.
  * (33) D. Peter, D. Komatitsch, Y. Luo, R. Martin, N. Le Goff, E. Casarotti, P. Le Loher, F. Magnoni, Q. Liu, C. Blitz, T. Nisson-Meyer, P. Basini, and J. Tromp, Forward and adjoint simulations of seismic wave propagation on fully unstructured hexahedral meshes, Geophysical Journal International, 186 (2011), pp. 721–739.
  * (34) J. Schütz and G. May, An adjoint consistency analysis for a class of hybrid mixed methods, IMA Journal of Numerical Analysis, (2013).
  * (35) C.-H. Teng, B.-Y. Lin, H.-C. Chang, H.-C. Hsu, C.-N. Lin, and K.-A. Feng, A Legendre pseudospectral penalty scheme for solving time-domain Maxwell’s equations, Journal of Scientific Computing, 36 (2008), pp. 351–390.
  * (36) F. Tröltzsch, Optimal Control of Partial Differential Equations: Theory, Methods and Applications, vol. 112 of Graduate Studies in Mathematics, American Mathematical Society, 2010.
  * (37) J. Utke, L. Hascoët, P. Heimbach, C. Hill, P. Hovland, and U. Naumann, Toward adjoinable MPI, in Parallel & Distributed Processing, 2009. IPDPS 2009. IEEE International Symposium on, IEEE, 2009, pp. 1–8.
  * (38) A. Walther, Automatic differentiation of explicit Runge-Kutta methods for optimal control, Computational Optimization and Applications, 36 (2007), pp. 83–108.  10.1007/s10589-006-0397-3.
  * (39) L. C. Wilcox, G. Stadler, C. Burstedde, and O. Ghattas, A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media, Journal of Computational Physics, 229 (2010), pp. 9373–9396.
  * (40) I. Yousept, Optimal control of Maxwell’s equations with regularized state constraints, Computational Optimization and Applications, 52 (2012), pp. 559–581.