# Ordinary differential equations for the adjoint Euler equations

Jacques Peter<EMAIL_ADDRESS>Jean-Antoine Désidéri  jean-
<EMAIL_ADDRESS>
###### Abstract

Ordinary Differential Equations are derived for the adjoint Euler equations
firstly using the method of characteristics in 2D. For this system of partial-
differential equations, the characteristic curves appear to be the
streamtraces and the well-known $\cal{C}^{+}$and $\cal{C}^{-}$curves of the
theory applied to the flow. The differential equations satisfied along the
streamtraces in 2D are then extended and demonstrated in 3D by linear
combinations of the original adjoint equations. These findings extend their
well-known counterparts for the direct system, and should serve analytical and
possibly numerical studies of the perfect-flow model with respect to adjoint
fields or sensitivity questions. Beside the analytical theory, the results are
demonstrated by the numerical integration of the compatibility relationships
for discrete 2D flow-fields and dual-consistent adjoint fields over a very
fine grid about an airfoil.

###### keywords:

continuous adjoint method, compressible Euler equations, supersonic flow,
charateristic curves

##  1  Introduction

In 1988, Jameson derived the continuous adjoint equations associated with the
2D and 3D Euler equations using general curvilinear coordinates [1]. With this
landmark article, the fluid dynamics and aeronautical communities became
better aware of the potential of the adjoint approach for design, that is, the
possibility to calculate gradient information at a cost scaling with the
number of functions to be differentiated, independently of the number of
design parameters. The equations in [1] appeared to be a natural starting
point for local optimizations involving a large number of design variables by
using adjoint gradients. However, in that setting, the flow and the dual
fields had to be calculated over a structured mesh.

Nine years later, Anderson and Venkatakrishnan [2, 3] and also Giles and
Pierce [4] derived the corresponding equations in Cartesian coordinates thus
allowing the application of the continuous approach (sometimes referred to as
the differentiate-then-discretize approach) on all types of meshes and, in
particular, on unstructured meshes. For the sake of simplicity, we present
here the two-dimensional case only in which the adjoint equations read

$-A^{T}\frac{\partial\psi}{\partial x}-B^{T}\frac{\partial\psi}{\partial
y}=0,\quad\mbox{in }\Omega\mbox{ the fluid domain}$ (1)

where $A$ and $B$ are the Jacobian matrices of the flux vectors $F_{x}$ and
$F_{y}$ of the Euler equations in the $x$ and $y$ directions respectively:

$\displaystyle F_{x}=\left(\begin{array}[]{c}\rho u\\\ \rho u^{2}+p\\\ \rho
uv\\\ \rho uH\\\ \end{array}\right)\qquad F_{y}=\left(\begin{array}[]{c}\rho
v\\\ \rho uv\\\ \rho v^{2}+p\\\ \rho vH\\\ \end{array}\right),$

with $\rho$ the density, $(u,v)$ the velocity components, $p$ the static
pressure and $H$ the total enthalpy. In the most common case where the
quantity of interest (QoI) is a line integral along the solid wall
$\Gamma_{w}$, it can be shown easily that the adjoint wall boundary condition
is well-posed provided that the function of interest depends only on the
static pressure [3]. In the classical case where the functional output of
interest is the force on $\Gamma_{w}$ projected in direction $\overline{d}$,
$J=\displaystyle\int_{\Gamma_{w}}p(\overline{n}\cdot\overline{d})ds$, the wall
boundary condition reads

$\overline{n}\cdot\overline{d}+\psi_{2}n_{x}+\psi_{3}n_{y}=0~{}~{}~{}~{}\textrm{on}~{}~{}\Gamma_{w}.$
(3)

For the farfield of an external flow, as well as for the inlet and outlet of
an internal flow, the boundary conditions are derived from the theory of local
one-dimensional characteristic decomposition [2, 5]. Here, the continuous
adjoint Euler equations and the associated boundary conditions are abbreviated
as (AE). Along with the growing use of the adjoint method for shape
optimization, goal oriented mesh adaptation and also meta-modelling, stability
or control, great effort is being devoted to gain understanding in the
mathematical properties of the (AE) solutions. The main results are summarized
here before discussing the characteristic relations for the (AE) system.

After the derivation of the (AE) equations, the first demonstrated property
was also due to Giles and Pierce [4]: in the common case where the function of
interest is an integral along the wall, the authors proved that the first and
last components of the adjoint vector $\psi$, associated with mass and energy
conservation, satisfy $\psi_{1}=H\psi_{4}$.

Besides, the integration by parts yielding (1) is not valid in the entire
fluid domain in the presence of flow discontinuities. After a series of works
dealing with the quasi-1D Euler equations – see [5, 6] and references therein
– Baeza et al. presented the equations complementing (1) along a shock line
[6] (denoted here $\Sigma$ as in the original reference). The new equations
are derived by introducing a complementary set of Lagrange multipliers,
multiplying the Rankine-Hugoniot conditions, viewed as constraints on
$\Sigma$. Finally, the continuity of the adjoint field $\psi$ along $\Sigma$
is established, although $\nabla\psi$ may be discontinuous across $\Sigma$, as
well as $\psi$ over $\Gamma_{w}\cap\Sigma$, and a so-called internal boundary
condition is derived:

$\displaystyle(\partial{\bf\psi}/\partial
t).(F_{x}t_{x}+F_{y}t_{y})=0\qquad\textrm{on}\qquad\Sigma$ (4)
$\displaystyle\textrm{with}~{}t~{}\textrm{the unit vector tangent to }\Sigma$

Coquel et al., Lozano and Renac [7, 8, 9] have derived additional
relationships by using (4), the jump operator applied to (1) across $\Sigma$
and the Rankine-Hugoniot equations.

The fact that $\psi_{1}=H\psi_{4}$ can be proven simply by forming the linear
combination of the first three lines of system (1) with coefficients
$(1,u/2,v/2)$. This yields
$\overline{U}.\nabla\psi_{1}-H\overline{U}.\nabla\psi_{4}=0$ (with
$\overline{U}=(u,v)$ the velocity vector). Note that this was also derived in
[4] by an approach based on physical source terms, constituting an important
analysis technique for the adjoint field of usual QoIs. In particular, this
method proved to be very fruitful to identify the zones where numerical
divergence of the adjoint vector is observed and mathematical divergence of
the solutions of (AE) is suspected. For the sake of clarity and brevity, we
restrict the present discussion to 2D flows about lifting airfoils, and to two
of these zones, namely the stagnation streamline and the wall, and to the lift
and drag as functions of interest.

More precisely, Giles and Pierce [4] introduced four physical punctual source
terms (or Green’s functions in the classical mathematical vocabulary) denoted
here $\delta R^{1}$, $\delta R^{2}$, $\delta R^{3}$, $\delta R^{4}$. These
terms are added to the right hand-side of the linearised Euler equations and
correspond respectively to (i) a mass source at fixed stagnation pressure
$p_{0}$ and enthalpy $H$ ; (ii) a normal force ; (iii) a change in $H$ at
fixed $p$ and $p_{0}$; and (iv) a change in $p_{0}$ at fixed $p$ and $H$. They
are linearly independent. (We refer to the original reference for the detailed
expression of these source terms.) The resulting changes in the QoI $J$,
$\delta J^{l}$, can be expressed as the integral over the domain of
$\psi\delta R^{l}$ that is, the value at the source location since $\delta
R^{l}$ is a Green’s function. These source terms also admit a physical
interpretation and their influence on the flow can be understood in terms of
mechanical principles, and sometimes even quantified finally providing insight
in the adjoint field [4].

It has been observed that the lift adjoint exhibits numerical divergence at
the stagnation streamline and at the wall at subcritical flow conditions. Also
the drag and lift adjoint of a transonic airfoil exhibit numerical divergence
at the same locations if the foot of at least one shock wave is located
strictly upwind the trailing edge – see [10, 11] and references therein.
Reference [9] includes a careful verification of this physical perturbations
approach applied to the discrete adjoint with a preliminary assessment of the
consistency between the linear (discrete adjoint) and the non-linear (flow
perturbation) evaluations of the $\delta J^{l}$. After this verification step,
the non-linear perturbed flow approach has been used (considering the physical
source terms point of view prior to the classical adjoint) and it appeared
that: (a) $\delta R^{4}$ is the only source term causing a numerical
divergence of $\delta J$ in the vicinity of the wall and stagnation streamline
; (b) in transonic condition, the numerical divergence of $\delta CLp^{4}$ and
$\delta CDp^{4}$ in these zones is mainly due to the displacement of the
shockfoot (or shock-feet if two shocks are not based at the trailing edge) ;
(c) this numerical divergence is transferred to the adjoint components via the
inverse matrix of the source terms ; (d) this does not necessarily prevent the
numerical satisfaction of the adjoint lift- (resp. drag-) boundary condition
at the wall (3) as the calculation of $\psi_{2}n_{x}+\psi_{3}n_{y}$ in this
approach involves the product of $\delta CLp^{4}$ (resp. $\delta CDp^{4}$) by
$(un_{x}+vn_{y})$.

The method of characteristics for 2D inviscid supersonic flow is a classical
method for deriving ordinary differential equations and, potentially, explicit
algebraic relations satisfied along two families of curves, denoted
$\cal{C}^{+}$(left running with respect to (w.r.t.) a streamline) and
$\cal{C}^{-}$(right running). Here, we recall the derivation of the continuous
equations and study their counterparts for the (AE) equations.
For the flow, the method starting point is the Cauchy problem posed for the
Euler equations: knowing the state variables along a fixed curve (L), is it
possible to calculate their partial derivatives, in both space directions, for
all points of (L) (this being a necessary condition for the flow calculation
in the fluid domain) ? The flow variables in two neighboring points of (L),
denoted here $a$ and $b$, are linked by fluid dynamics equations and basic
first order Taylor formulas. The more general presentations deal with
rotational flows and, geometrically, both planar and axisymmetrical flows [12,
13, 14, 15]. The authors derive a system of equations for the derivatives of
selected primitive variables. This system is linear in the unknown partial
derivatives with non linear functions of the state variables as coefficients.
When its determinant is equal to zero, it cannot be solved. That is the case
if $b$ is on the streamline of $a$ or, in case of a supersonic flow, if the
angle of $\overrightarrow{ab}=(dx,dy)$ w.r.t. the streamline passing through
$a$ is $\pm\sin^{-1}(1/M)$ ($M$ being the local Mach number). The specific
curves (L) where these conditions are satisfied for every points are, for all
Mach numbers, the streamtraces, and, where the flow is supersonic, the so-
called $\cal{C}^{+}$(left running curves w.r.t. the streamtraces with angle
$\sin^{-1}(1/M)$) and the $\cal{C}^{-}$(right running curves the w.r.t.
streamtraces with angle $-\sin^{-1}(1/M)$). Along these curves, the physical
existence and boundedness of the vector of unknowns allows to conclude from
the nullity of the determinant in the denominator of Cramer’s formulas, to the
nullity of the determinants appearing in the numerators, and this, for all the
variables. The classical computational method for supersonic flow [12, 13, 16,
14, 17, 15] is supported by the corresponding differential forms valid along
the $\cal{C}^{+}$and $\cal{C}^{-}$curves, and the property of constant total
enthalpy and constant entropy along the streamtraces. If the flow is
irrotational, homoenthalpic and homoentropic, simpler equations are derived
for the velocity magnitude and the velocity angle [13, 18, 16] or the velocity
potential [17] and the differential forms satisfied along the $\cal{C}^{+}$and
$\cal{C}^{-}$curves may be integrated. This permits to establish the well-
known equations

$k^{-}=\phi+\nu(M)~{}~{}\textrm{is constant along a}~{}~{}{\cal
C}^{-}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}k^{+}=\phi-\nu(M)~{}~{}\textrm{is
constant along a}~{}~{}{\cal C}^{+}$ (5)

in which $\phi$ is the streamline angle, $\nu(M)$ the Prandtl-Meyer function,

$\phi=\tan^{-1}(v/u),\quad\nu(M)=\sqrt{\frac{\gamma+1}{\gamma-1}}\tan^{-1}(\sqrt{\frac{\gamma-1}{\gamma+1}(M^{2}-1)})-\tan^{-1}(\sqrt{M^{2}-1}),$

and $\gamma$ the ratio of specific heats ($\gamma=\frac{7}{5}$ for diatomic
perfect gas). Finally, let us recall that the integration of the corresponding
ODEs along the trajectories results in the property of constant enthalpy and
constant entropy.
Besides, Bonnet and Luneau indicate that the mechanical equations posed in $a$
may be expressed in the Cartesian frame of reference rather than in the usual
local frame derived from the local velocity [14]. Note that the 2D
characteristic equations could also be calculated in the Cartesian frame and
in an inexpert way, without taking advantage of the known properties of the
streamtraces. Then the following the $8\times 8$ linear system relating the
derivatives of the conservative variables would be solved:

$\begin{bmatrix}dx&0&0&0&dy&0&0&0\\\ 0&dx&0&0&0&dy&0&0\\\ 0&0&dx&0&0&0&dy&0\\\
0&0&0&dx&0&0&0&dy\\\ &&&&&&&\\\ &&&&&&&\\\ &&A&&&&B&\\\
&&&&&&&\end{bmatrix}\begin{bmatrix}(\partial\rho/\partial x)\\\ (\partial\rho
u/\partial x)\\\ (\partial\rho v/\partial x)\\\ (\partial\rho E/\partial x)\\\
(\partial\rho/\partial y)\\\ (\partial\rho u/\partial y)\\\ (\partial\rho
v/\partial y)\\\ (\partial\rho E/\partial
y)\end{bmatrix}=\begin{bmatrix}\rho^{b}-\rho^{a}\\\ \rho u^{b}-\rho u^{a}\\\
\rho v^{b}-\rho v^{a}\\\ \rho E^{b}-\rho E^{a}\\\ 0.\\\ 0.\\\ 0.\\\
0.\end{bmatrix}$ (6)

The starting point of our analytical development resides in the observation
that (6) and the corresponding linear system for the (AE) equations, (7), have
the same determinant. From this observation, the adjoint Euler characteristics
equations are established in Sec. 2. The theoretical findings are linked with
former researches and illustrated by numerical computational solutions over a
very fine grid in Sec. 3. Conclusions are drawn in Sec. 4.

##  2  Adjoint characteristic equations for 2D supersonic flow

The method exposed in [17] (resp. [15]) for potential (resp. general) inviscid
flow has served as a guideline to our derivation for the adjoint system. For
all our calculations, we assume an ideal gas law for the static pressure
$p=(\gamma-1)\rho e=(\gamma-1)(\rho E-0.5\rho||\overline{U}||^{2})$ with a
constant $\gamma$.

###  2.1 Problem statement

Given two fixed close points in the supersonic zone, $a$ and $b$, is it
possible to estimate $(\partial\psi/\partial x),(\partial\psi/\partial y)$
from the local value of the flow field and ($\psi_{a}$,$\psi_{b}$) ? This
question is the starting point of the method of characteristics in which
specific lines are identified along which this problem is ill-posed, and
particular ordinary differential equations are satisfied. Let us denote
$\overrightarrow{ab}=(dx,dy)$ and first assume that $dx\neq 0$. By definition
of differential forms, and in view of the adjoint system (1), the following
holds

$\begin{bmatrix}dx&0&0&0&dy&0&0&0\\\ 0&dx&0&0&0&dy&0&0\\\ 0&0&dx&0&0&0&dy&0\\\
0&0&0&dx&0&0&0&dy\\\ &&&&&&&\\\ &&&&&&&\\\ &-&A^{T}&&&-&B^{T}&\\\
&&&&&&&\end{bmatrix}\begin{bmatrix}(\partial\psi_{1}/\partial x)\\\
(\partial\psi_{2}/\partial x)\\\ (\partial\psi_{3}/\partial x)\\\
(\partial\psi_{4}/\partial x)\\\ (\partial\psi_{1}/\partial y)\\\
(\partial\psi_{2}/\partial y)\\\ (\partial\psi_{3}/\partial y)\\\
(\partial\psi_{4}/\partial y)\end{bmatrix}=\begin{bmatrix}d\psi_{1}\\\
d\psi_{2}\\\ d\psi_{3}\\\ d\psi_{4}\\\ 0\\\ 0\\\ 0\\\ 0\end{bmatrix}$ (7)

in which by neglecting second-order terms in space:
$(d\psi_{1},d\psi_{2},d\psi_{3},d\psi_{4})=(\psi^{b}_{1}-\psi^{a}_{1},\psi^{b}_{2}-\psi^{a}_{2},\psi^{b}_{3}-\psi^{a}_{3},\psi^{b}_{4}-\psi^{a}_{4})$.
The determinant of the linear system is evidently

$\begin{vmatrix}&&&&&&\\\ &dxI&&&&dyI&\\\ &&&&&&\\\ &&&&&&\\\ &&&&&&\\\
&-A^{T}&&&&-B^{T}&\\\ &&&&&&\end{vmatrix}=\begin{vmatrix}&&&&&&\\\
&dxI&&&0&&\\\ &&&&&&\\\ &&&&&&\\\ &&&&&&\\\ &-A^{T}&&&-B^{T}+dy/dxA^{T}&&\\\
&&&&&&\end{vmatrix}=dx^{4}|-B^{T}+dy/dxA^{T}|=|-dxB^{T}+dyA^{T}|$

Of course, $|-dxB^{T}+dyA^{T}|$ is equal to $|-dxB+dyA|$ and the value of this
determinant is known from the eigenvalues of the matrix:

$D=|-dxB+dyA|=(-v~{}dx+u~{}dy)^{2}(-v~{}dx+u~{}dy+c~{}ds)(-v~{}dx+u~{}dy-c~{}ds),$

in which

$c=\sqrt{\frac{\gamma p}{\rho}},\quad ds=\sqrt{dx^{2}+dy^{2}}.$

Similarly to the flow derivatives reconstruction [17, 15], the problem of
adjoint derivatives reconstruction in a supersonic zone is ill-posed along the
same three families of curves

$\displaystyle-v~{}dx+u~{}dy$ $\displaystyle=$ $\displaystyle
0\qquad\qquad\textrm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$\cal{S}$~{}
streamtraces \qquad(all Mach numbers)}}$ (8)
$\displaystyle-v~{}dx+u~{}dy+c~{}ds$ $\displaystyle=$ $\displaystyle
0\qquad\qquad\textrm{$\cal{C}^{-}$characteristics
\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(supersonic
flow only)}}$ (9) $\displaystyle-v~{}dx+u~{}dy-c~{}ds$ $\displaystyle=$
$\displaystyle 0\qquad\qquad\textrm{$\cal{C}^{+}$characteristics
\quad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(supersonic
flow only)}}$ (10)

Figure 1: $\cal{S}$, $\cal{C}^{+}$and $\cal{C}^{-}$curves according to local
Mach number. Neighboring $a$ and $b$ points along these curves

Classically, the method of characteristics uses the ill-posedness of (7) in
the following way: along the curves defined by equations (8), (9) or (10), not
only the denominator appearing in the Cramer formulas applied to the linear
equations (7) is equal to zero, but the numerators giving the eight components
of $(\partial\psi/\partial x)$ and $(\partial\psi/\partial y)$ must also be
equal to zero for the fractions not to be singular. This (somehow paradoxical)
technique allows the derivation of equations (5). It is derived here for the
(AE) system by analysing the set of linear equations (7).

###  2.2 Null differential forms in the adjoint variations along trajectories
and characteristics

The transposed of the Euler flux Jacobian matrices in $x$ and $y$ direction
read

$A^{T}=\begin{bmatrix}0&\gamma_{1}E_{c}-u^{2}&-uv&(\gamma_{1}E_{c}-H)u\\\
1&(3-\gamma)u&v&H-\gamma_{1}u^{2}\\\ 0&-\gamma_{1}v&u&-\gamma_{1}uv\\\
0&\gamma_{1}&0&\gamma u\\\
\end{bmatrix}~{}~{}~{}~{}~{}~{}~{}~{}B^{T}=\begin{bmatrix}0&-uv&\gamma_{1}E_{c}-v^{2}&(\gamma_{1}E_{c}-H)v\\\
0&v&-\gamma_{1}u&-\gamma_{1}uv\\\ 1&u&(3-\gamma)v&H-\gamma_{1}v^{2}\\\
0&0&\gamma_{1}&\gamma v\\\ \end{bmatrix}$

in the usual notations in aerodynamics and $\gamma_{1}=\gamma-1$. Let

$t=\frac{dy}{dx},\quad\kappa=ut-v,$

and also introduce the following notations for the column vectors of the
transposed Jacobian matrices: $A^{T}=[A_{1}|A_{2}|A_{3}|A_{4}]$,
$B^{T}=[B_{1}|B_{2}|B_{3}|B_{4}]$. Before presenting the results, the
principle of the calculation is recalled in one of the cases that leads to the
simplest calculations: the definition of $(\partial\psi_{4}/\partial x)$ along
the curves where (8), (9) or (10) is satisfied (that is, the streamtraces, the
$\cal{C}^{-}$or $\cal{C}^{+}$characteristic) requires that, along these curves

$\begin{vmatrix}~{}dx&0&0&d\psi_{1}&dy&0&0&0~{}\\\
~{}0&dx&0&d\psi_{2}&0&dy&0&0~{}\\\ ~{}0&0&dx&d\psi_{3}&0&0&dy&0~{}\\\
~{}0&0&0&d\psi_{4}&0&0&0&dy~{}\\\ ~{}|&|&|&|&|&|&|&|~{}\\\
~{}-A_{1}&-A_{2}&-A_{3}&0&-B_{1}&-B_{2}&-B_{3}&-B_{4}~{}\\\
~{}|&|&|&|&|&|&|&|~{}\\\ \end{vmatrix}=0$

The determinant is expanded along the fourth column and the following
notations are used

$-C_{4x}^{1}d\psi_{1}+C_{4x}^{2}d\psi_{2}-C_{4x}^{3}d\psi_{3}+C_{4x}^{4}d\psi_{4}=0$
(11)

in which, for example

$C_{4x}^{1}=\begin{vmatrix}~{}0&dx&0&0&dy&0&0~{}\\\ ~{}0&0&dx&0&0&dy&0~{}\\\
~{}0&0&0&0&0&0&dy~{}\\\ ~{}|&|&|&|&|&|&|~{}\\\
~{}-A_{1}&-A_{2}&-A_{3}&-B_{1}&-B_{2}&-B_{3}&-B_{4}~{}\\\
~{}|&|&|&|&|&|&|~{}\\\ \end{vmatrix}=\begin{vmatrix}~{}0&dx&0&0&0&0&0~{}\\\
~{}0&0&dx&0&0&0&0~{}\\\ ~{}0&0&0&0&0&0&dy~{}\\\ ~{}|&|&|&|&|&|&|~{}\\\
~{}-A_{1}&-A_{2}&-A_{3}&-B_{1}&-B_{2}+tA_{2}&-B_{3}+tA_{3}&-B_{4}~{}\\\
~{}|&|&|&|&|&|&|~{}\\\ \end{vmatrix}$

Finally

$C_{4x}^{1}=dx^{2}dy~{}|~{}~{}-A_{1}~{}~{}-B_{1}~{}~{}(-B_{2}+tA_{2})~{}~{}(-B_{3}+tA_{3})~{}~{}|$

The determinant of this $4\times 4$ matrix is easily calculated thanks to the
simplicity of the first two columns $A_{1}$ and $B_{1}$. The final result is

$C_{4x}^{1}=dx^{2}dy~{}\gamma_{1}~{}\kappa~{}(u+vt)$

We emphasise that at this stage no assumption is made on the value of
$t=dy/dx$ w.r.t. the velocity vector $(u,v)$. In particular $t$ is not assumed
to be the tangent of the angle of the velocity w.r.t. the $x$ axis and
$\kappa$ is not assumed to be zero. This is mandatory to derive relations that
can be used for all three types of specific curves and also to account for the
multiplicity of the eigenvalue $(-v~{}dx+u~{}dy)$ along streamtraces. The
other terms of the differential form of interest read

$\displaystyle C_{4x}^{2}$ $\displaystyle=$ $\displaystyle-
dx^{2}dy~{}|~{}~{}-A_{2}~{}~{}(-B_{1}+tA_{1})~{}~{}-B_{2}~{}~{}(-B_{3}+tA_{3})~{}~{}|=-dx^{2}dy~{}\gamma_{1}~{}\kappa~{}(u^{2}+v^{2})$
$\displaystyle C_{4x}^{3}$ $\displaystyle=$ $\displaystyle
dx^{2}dy~{}|~{}~{}-A_{3}~{}~{}(-B_{1}+tA_{1})~{}~{}(-B_{2}+tA_{2})~{}~{}-B_{3}~{}~{}|=dx^{2}dy~{}\gamma_{1}~{}\kappa~{}t~{}(u^{2}+v^{2})$
$\displaystyle C_{4x}^{4}$ $\displaystyle=$ $\displaystyle
dx^{3}~{}|~{}~{}(-B_{1}+tA_{1})~{}~{}(-B_{2}+tA_{2})~{}~{}(-B_{3}+tA_{3})~{}~{}-B_{4}~{}~{}|$
$\displaystyle=$ $\displaystyle-dx^{3}\kappa\left((\gamma_{1}+\gamma
t^{2})u^{2}v-2uv^{2}t+\gamma_{1}H\kappa+(\gamma+\gamma_{1}t^{2})v^{3}-\gamma_{1}v~{}(1+t^{2})~{}E_{c}\right)$

The explicit expression of (11), the necessary condition for the boundedness
of $(\partial\psi_{4}/\partial x)$, hence reads

$\displaystyle-$ $\displaystyle
dx^{2}dy~{}\gamma_{1}~{}\kappa~{}(u+vt)~{}d\psi_{1}-dx^{2}dy~{}\gamma_{1}~{}\kappa~{}(u^{2}+v^{2})d\psi_{2}-dx^{2}dy~{}\gamma_{1}~{}\kappa~{}t~{}(u^{2}+v^{2})d\psi_{3}$
(12) $\displaystyle-$ $\displaystyle dx^{3}\kappa\left((\gamma_{1}+\gamma
t^{2})u^{2}v-2uv^{2}t+\gamma_{1}H\kappa+(\gamma+\gamma_{1}t^{2})v^{3}-\gamma_{1}v~{}E_{c}~{}(1+t^{2})\right)d\psi_{4}=0$

Assuming that $dx\neq 0$, this equation may be further simplified for the
$\cal{C}^{-}$and $\cal{C}^{+}$for which $\kappa\neq 0$:

$\displaystyle\gamma_{1}~{}t(u+vt)~{}d\psi_{1}+\gamma_{1}~{}t(u^{2}+v^{2})~{}d\psi_{2}+~{}\gamma_{1}~{}t^{2}~{}(u^{2}+v^{2})~{}d\psi_{3}$
(13) $\displaystyle+$ $\displaystyle\left((\gamma_{1}+\gamma
t^{2})u^{2}v-2uv^{2}t+\gamma_{1}H\kappa+(\gamma+\gamma_{1}t^{2})v^{3}-\gamma_{1}v~{}E_{c}~{}(1+t^{2})\right)d\psi_{4}=0.$

As $(-v~{}dx+u~{}dy)=\kappa~{}dx$ has a multiplicity of two in the determinant
of (7), equation (13) is also needed for the existence of
$(\partial\psi_{4}/\partial x)$ and hence true for neighboring points $a$ and
$b$ along the same $\cal{S}$curves. (This point is detailed in §2.3.)
For the sake of clarity, the results of the corresponding calculations for the
existence of the seven other partial derivatives along the $\cal{S}$,
$\cal{C}^{-}$and $\cal{C}^{+}$curves, and the properties of the $C^{i}_{jx}$
$C^{i}_{jy}$ coefficients are presented in Appendix A and B. Only the
counterparts of equation (13), for the existence of
$(\partial\psi_{1}/\partial x)$, $(\partial\psi_{2}/\partial x)$,
$(\partial\psi_{3}/\partial x)$, $(\partial\psi_{1}/\partial y)$,
$(\partial\psi_{2}/\partial y)$, $(\partial\psi_{3}/\partial y)$ and
$(\partial\psi_{4}/\partial y)$ are presented hereafter in this order:

$\displaystyle((2tu+(t^{2}-1)v)(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)-(u+vt)(\gamma_{1}tH+\gamma_{1}u\kappa+\gamma v\kappa t)d\psi_{1}$
(14) $\displaystyle+$
$\displaystyle~{}\gamma_{1}~{}t~{}(u^{2}+v^{2})~{}H(d\psi_{2}+td\psi_{3})+\gamma_{1}~{}t~{}(u+vt)~{}H^{2}d\psi_{4}=0$
$\displaystyle t~{}(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)~{}(d\psi_{1}+Hd\psi_{4})$ (15) $\displaystyle-$
$\displaystyle(\gamma_{1}u^{2}v-uv^{2}t+\gamma
v^{3}-\gamma_{1}(ut+v)(E_{c}+H))~{}(d\psi_{2}+td\psi_{3})=0$ $\displaystyle
t~{}(t\gamma_{1}H+t\gamma_{1}E_{c}-\gamma
u\kappa)(d\psi_{1}+Hd\psi_{4})+t~{}((\gamma+1)u^{2}v-t\gamma_{1}uv^{2}+(v+ut)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2}))d\psi_{2}$ (16) $\displaystyle+$
$\displaystyle~{}((tv^{2}-u\kappa)(-\gamma_{1}u+(\gamma+1)tv)-(ut-(1+2t^{2})v)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
v^{2}))d\psi_{3}=0$ $\displaystyle\left((\gamma_{1}+\gamma
t^{2})u^{3}-2u^{2}vt-\gamma_{1}t\kappa
H+(\gamma+\gamma_{1}t^{2})uv^{2}-\gamma_{1}u~{}(1+t^{2})~{}E_{c}\right)d\psi_{1}$
(17) $\displaystyle+$
$\displaystyle~{}\gamma_{1}~{}(u^{2}+v^{2})~{}H(d\psi_{2}+td\psi_{3})+~{}\gamma_{1}~{}(u+vt)~{}H^{2}d\psi_{4}=0$
$\displaystyle(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)~{}(d\psi_{1}+Hd\psi_{4})$ (18) $\displaystyle+$
$\displaystyle((v\kappa+u^{2})(vt-(\gamma_{1}+\gamma
t^{2})u)-(vt-(2+t^{2})u)(\gamma_{1}E_{c}+\gamma_{1}H+\gamma
v\kappa))d\psi_{2}$ $\displaystyle-$
$\displaystyle(\gamma_{1}u^{2}v-uv^{2}t+\gamma
v^{3}-\gamma_{1}(ut+v)(E_{c}+H))d\psi_{3}=0$
$\displaystyle~{}(t\gamma_{1}H+t\gamma_{1}E_{c}-\gamma
u\kappa)(d\psi_{1}+Hd\psi_{4})$ (19) $\displaystyle+$
$\displaystyle((\gamma+1)u^{2}v-t\gamma_{1}uv^{2}+(v+ut)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2}))(d\psi_{2}+td\psi_{3})=0$
$\displaystyle~{}\gamma_{1}~{}(u+vt)~{}d\psi_{1}+\gamma_{1}~{}(u^{2}+v^{2})(d\psi_{2}+td\psi_{3})$
(20) $\displaystyle+$ $\displaystyle\left((\gamma_{1}+\gamma
t^{2})u^{3}-2u^{2}vt-\gamma_{1}t\kappa
H+(\gamma+\gamma_{1}t^{2})uv^{2}-\gamma_{1}u~{}E_{c}~{}(1+t^{2})\right)d\psi_{4}=0$

###  2.3 Ordinary differential equations for the adjoint along the
streamtraces $\cal{S}$

The trajectories are one of the families of specific curves for the gradient
calculation problem (7). Along these curves $u~{}dy-v~{}dx=0$ is a zero of the
denominator of Cramer’s formulas with multiplicity two. Let us first assume
that point $a$ is fixed and point $b$ is very close to ${\cal S}_{a}$, the
streamtrace passing through $a$ but not on this curve. The first-order
expression of $(\partial\psi_{1}/\partial x)$ is

$\frac{\partial\psi_{1}}{\partial
x}=\frac{C^{1}_{1x}d\psi_{1}-C^{2}_{1x}d\psi_{2}+C^{3}_{1x}d\psi_{3}-C^{4}_{1x}d\psi_{4}}{(-v~{}dx+u~{}dy)^{2}(-v~{}dx+u~{}dy+c~{}ds)(-v~{}dx+u~{}dy-c~{}ds)}$

Actually $\kappa dx=u~{}dy-v~{}dx$ is a factor of all four coefficients
$C^{1}_{1x}$, $C^{2}_{1x}$, $C^{3}_{1x}$ and $C^{4}_{1x}$. We denote by
$\bar{C}^{l}_{mx}$ the coefficients obtained by removing the $(\kappa~{}dx)$
factor from the corresponding $C^{l}_{mx}$. Obviously

$\frac{\partial\psi_{1}}{\partial
x}=\frac{\bar{C}^{1}_{1x}d\psi_{1}-\bar{C}^{2}_{1x}d\psi_{2}+\bar{C}^{3}_{1x}d\psi_{3}-\bar{C}^{4}_{1x}d\psi_{4}}{(-v~{}dx+u~{}dy)(-v~{}dx+u~{}dy+c~{}ds)(-v~{}dx+u~{}dy-c~{}ds)}$

If point $b$ is moved closer and closer to ${\cal S}_{a}$
$(-v~{}dx+u~{}dy)\rightarrow 0$, so that the boundedness of
$(\partial\psi_{1}/\partial x)$ requires that

$\bar{C}^{1}_{1x}d\psi_{1}-\bar{C}^{2}_{1x}d\psi_{2}+\bar{C}^{3}_{1x}d\psi_{3}-\bar{C}^{4}_{1x}d\psi_{4}=0~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{on}~{}~{}~{}{\cal
S}_{a}$

This expression and its counterparts for the other derivatives
$(\partial\psi_{2}/\partial x)$ … $(\partial\psi_{4}/\partial y)$ have to be
satisfied for all trajectories. How many of these eight differential forms are
independent ? If $\kappa~{}dx=0$, then

$\bar{C}^{1}_{1x}=-t~{}\bar{C}^{1}_{1y}~{}~{}~{}~{}~{}\bar{C}^{2}_{2x}=-t~{}\bar{C}^{2}_{2y}~{}~{}~{}~{}\bar{C}^{3}_{3x}=-t~{}\bar{C}^{3}_{3y}~{}~{}~{}~{}\bar{C}^{4}_{4x}=-t~{}\bar{C}^{4}_{4y}$
(21)

as in this case

$~{}\bar{C}^{1}_{1x}=\bar{C}^{4}_{4x}=-0.5~{}t\gamma_{1}dx^{2}(1+t^{2})^{2}u^{3}~{}~{}~{}~{}~{}~{}~{}\bar{C}^{1}_{1y}=\bar{C}^{4}_{4y}=0.5~{}\gamma_{1}dx^{2}(1+t^{2})^{2}u^{3}$
$~{}\bar{C}^{2}_{2x}=2~{}dx^{2}\gamma_{1}tuH~{}~{}~{}~{}~{}\bar{C}^{2}_{2y}=-2~{}dx^{2}\gamma_{1}uH~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bar{C}^{3}_{3x}=2~{}dx^{2}\gamma_{1}t^{3}uH~{}~{}~{}~{}~{}\bar{C}^{3}_{3y}=-2~{}dx^{2}\gamma_{1}t^{2}uH.$

Relations (40) to (43) are valid for the $\bar{C}$ coefficients (as they stand
whatever the values of $w$ and $dx$, they may be simplified by $wdx$). In the
specific case where $\kappa=0~{}$ they are completed by (21). Equations
{(17),(18),(19),(20)} (necessary for the boundedness of the
$\partial\psi^{l}/\partial y$) and {(14),(15),(16),(13)} (same for
$\partial\psi^{l}/\partial x$) are then proportional by a $(-t)$ factor.
Considering the range of the set of the eight differential forms, it appears
that one of these two sets of four equations need be accounted for.

The relations stemming from the existence of the $\partial x$ partial
derivative are retained. Equation (14), required for the definition of
$(\partial\psi_{1}/\partial x)$, is further simplified using the specific
properties of a trajectory ($\kappa=0~{}~{}v=ut$):

$0.5t\gamma_{1}(1+t^{2})^{2}u^{3}d\psi_{1}+~{}\gamma_{1}~{}t~{}(1+t^{2})u^{2}~{}H(d\psi_{2}+td\psi_{3})+~{}\gamma_{1}~{}t~{}(1+t^{2})u~{}H^{2}d\psi_{4}=0$
$0.5(1+t^{2})u^{2}d\psi_{1}+u~{}H(d\psi_{2}+td\psi_{3})+~{}H^{2}d\psi_{4}=0$

For the streamstraces, the equation finally derived from the existence of
$(\partial\psi_{1}/\partial x)$ is

$E_{c}~{}d\psi_{1}+H(u~{}d\psi_{2}+v~{}d\psi_{3})+~{}H^{2}d\psi_{4}=0$ (22)

Note that we have assumed that $u\neq 0$ and $dx\neq 0$ to perform the
calculations but finally obtained an expression that is also well-defined in
this specific case. Equation (15), is further simplified for the motion along
a trajectory:

$~{}t~{}(\gamma_{1}H+\gamma_{1}E_{c})~{}(d\psi_{1}+Hd\psi_{4})+2\gamma_{1}~{}u~{}t~{}H~{}(d\psi_{2}+td\psi_{3})=0$
$(H+E_{c})~{}(d\psi_{1}+Hd\psi_{4})+2H~{}(u~{}d\psi_{2}+v~{}d\psi_{3})=0$

Using the first relation and simplifying by $H$, we get

$H(d\psi_{1}+Hd\psi_{4})+2H~{}(u~{}d\psi_{2}+v~{}d\psi_{3})+E_{c}Hd\psi_{4}-H(ud\psi_{2}+vd\psi_{3})-~{}H^{2}d\psi_{4}=0$
$\textrm{and
finally}~{}~{}~{}~{}~{}~{}~{}d\psi_{1}+u~{}d\psi_{2}+v~{}d\psi_{3}+E_{c}d\psi_{4}=0$
(23)

Simplifying equation (16) for trajectories yields

$(H+E_{c})(d\psi_{1}+Hd\psi_{4})+2H(ud\psi_{2}+vd\psi_{3})=0,$

that had already been derived. Using $\kappa=0$ and $v=ut$, the fourth
relation (first simplified by a $\kappa dx$ factor) gives

$\gamma_{1}~{}tu(1+t^{2})~{}d\psi_{1}+\gamma_{1}~{}tu^{2}(1+t^{2})~{}d\psi_{2}+~{}\gamma_{1}t^{2}~{}u^{2}(1+t^{2})~{}d\psi_{3}+0.5t\gamma_{1}(1+t^{2})^{2}u^{3}d\psi_{4}=0$
$\textrm{or}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}d\psi_{1}+u~{}d\psi_{2}+tu~{}d\psi_{3}+0.5(1+t^{2})u^{2}d\psi_{4}=0$
$\textrm{and
finally}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}d\psi_{1}+u~{}d\psi_{2}+v~{}d\psi_{3}+E_{c}~{}d\psi_{4}=0$

This equation is similar to (23). These differential forms along the
trajectories ${\cal S}$ may be turned in differential equations. The natural
variable w.r.t. which differentiate, is the curvilinear abscissa along the
streamtraces, $s$, increasing in the direction opposite to the motion (as this
is the direction of the adjoint transport of information from the support of
the function of interest). The final equations along the ${\cal S}$ curves are
then

$\displaystyle
E_{c}~{}\frac{d\psi_{1}}{ds}+H(u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds})+~{}H^{2}\frac{d\psi_{4}}{ds}=0$
(24)
$\displaystyle\frac{d\psi_{1}}{ds}+u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+E_{c}\frac{d\psi_{4}}{ds}=0$
(25)

Although the calculations of this section are not very complex, they have been
checked with the computer algebra software Maple. The independent variables of
the formal calculations are $(u,v)$ also $M$ the Mach number (that allows the
calculation of the speed of sound $c$ and then the total enthalpy $H$) and
$\gamma$. The set of four differential forms {(17), (18), (19), (20)} was
associated with {(22),(23)} in a $6\times 4$ matrix that was again found to
have rank two.

###  2.4 Ordinary differential equation for the adjoint along the
$\cal{C}^{+}$and $\cal{C}^{-}$characteristics

Let us first note that the determinant in the denominator of the Cramer
formulas for (7) may be calculated as

$D=dx~{}C^{1}_{1x}+dy~{}C^{1}_{1y}$

developing $D$ along the first line. Doing the same along the second, third
and fourth lines yields

$D=dx~{}C^{2}_{2x}+dy~{}C^{2}_{2y}=dx~{}C^{3}_{3x}+dy~{}C^{3}_{3y}=dx~{}C^{4}_{4x}+dy~{}C^{4}_{4y}$

As $D=0$ along the $\cal{C}^{+}$and $\cal{C}^{-}$characteristics, equations
(40) to (43) are completed by

$C^{1}_{1x}=-t~{}C^{1}_{1y}\qquad C^{2}_{2x}=-t~{}C^{2}_{2y}\qquad
C^{3}_{3x}=-t~{}C^{3}_{3y}\qquad C^{4}_{4x}=-t~{}C^{4}_{4y}$

so that the differential forms stemming from the boundedness of
$(\partial\psi^{l}/\partial x)$ and their counterparts for
$(\partial\psi^{l}/\partial x)$ are proportional. (Incidentally, note that
this argument may have been used also in the previous subsection.)
Numerical tests indicate that the four differential forms associated with the
existence of (choosing the second set) $(\partial\psi_{1}/\partial y)$,
$(\partial\psi_{2}/\partial y)$, $(\partial\psi_{3}/\partial y)$,
$(\partial\psi_{4}/\partial y)$ are proportional but the corresponding
calculations are much more complex than in the previous subsection as the
expression of $t$ is now:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t^{\pm}}=\tan(\phi+\beta)~{}~{}~{}~{}\textrm{with}~{}~{}~{}\tan\phi=\frac{v}{u}~{}~{}~{}\sin\beta=\pm\frac{1}{M}$
(26)

with, of course $\sin\beta=1/M$ for the $\cal{C}^{+}$and $\sin\beta=-1/M$ for
the $\cal{C}^{-}$curves.
With this expression of $t$ and $\kappa$ not being equal to zero, searching
the rank of {(17), (18), (19), (20)} is much more difficult. However, the task
was successfully accomplished with the assistance of Maple, using once again
$(u,v,M,\gamma)$ as independent formal variables. Rank one was indicated by
Maple and the correspondence with the counterpart characteristics for the flow
seemed very sound. Knowing this result from formal calculation, its
demonstration by hand was searched for.
First, the value of $t$ is calculated along the $\cal{C}^{+}$and
$\cal{C}^{-}$curves. For $\beta$ in $[-\pi/2,\pi/2]$ (26) yields

$t^{\pm}=\frac{\displaystyle\frac{v}{u}\pm\displaystyle\frac{1}{\sqrt{M^{2}-1}}}{1\mp\displaystyle\frac{v}{u}\displaystyle\frac{1}{\sqrt{M^{2}-1}}}=\frac{uvM^{2}\pm(u^{2}+v^{2})\sqrt{M^{2}-1}}{u^{2}(M^{2}-1)-v^{2}}=\frac{uv\pm
c^{2}\sqrt{M^{2}-1}}{u^{2}-c^{2}}=\frac{uv\pm
c\sqrt{u^{2}+v^{2}-c^{2}}}{u^{2}-c^{2}}$ (27)

We then note that if the differential forms (17) and (20) were proportional,
from the expressions of $C^{2}_{1y},C^{3}_{1y},C^{2}_{4y}$ and $C^{3}_{4y}$,
the ratio between their terms would be $(-H)$. It is then easily checked that
the two non-trivial conditions for proportionality, $C^{1}_{1y}=-HC^{4}_{1y}$
and $C^{1}_{4y}=-HC^{4}_{4y}$, are equivalent to a single equality

$\gamma_{1}~{}(u+vt)H=\left((\gamma_{1}+\gamma
t^{2})u^{3}-2u^{2}vt-\gamma_{1}t\kappa
H+(\gamma+\gamma_{1}t^{2})uv^{2}-\gamma_{1}u~{}(1+t^{2})~{}E_{c}\right)$ (28)

Wherever $u\neq 0$, this condition is equivalent to

$\gamma_{1}~{}(1+t^{2})~{}H=\gamma_{1}~{}(1+t^{2})~{}E_{c}+(tu-v)^{2}$ (29)

that is, precisely, the degree two equation which roots are the values of $t$
along the $\cal{C}^{+}$and $\cal{C}^{-}$curves (27). Along these curves, using
(28) that is now an established property along the $\cal{C}^{+}$and
$\cal{C}^{-}$, these two differential forms may be simplified as

$(u+vt^{\pm})~{}d\psi_{1}+(u^{2}+v^{2})(d\psi_{2}+t^{\pm}d\psi_{3})+H(u+vt^{\pm})~{}d\psi_{4}=0,$
(30)

or, under the form of an ordinary differential equation,

$(u+vt^{\pm})~{}\frac{d\psi_{1}}{ds}+(u^{2}+v^{2})(\frac{d\psi_{2}}{ds}+t^{\pm}\frac{d\psi_{3}}{ds})+H(u+vt^{\pm})~{}\frac{d\psi_{4}}{ds}=0.$
(31)

Comparing equations (19) – expressing the boundedness of
$(\partial\psi_{3}/\partial y)$ – and (30) – its counterpart for
$(\partial\psi_{1}/\partial y)$ and $(\partial\psi_{4}/\partial y)$) – it is
easily derived that these equations are proportional on the $\cal{C}^{+}$and
$\cal{C}^{-}$curves if and only if for $t=t^{\pm}$

$(t\gamma_{1}H+t\gamma_{1}E_{c}-\gamma
u\kappa)(u^{2}+v^{2})=((\gamma+1)u^{2}v-t\gamma_{1}uv^{2}+(v+ut)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2}))(u+vt)$

Wherever $uv\neq 0$, this equation is equivalent to

$\gamma(u^{2}+v^{2})=((\gamma+1)u-t\gamma_{1}v)(u+vt)+(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2})(1+t^{2})$

that is found to be equivalent to equation (29), the degree two equation $t$
which roots are the slope coefficients of the $\cal{C}^{+}$and $\cal{C}^{-}$.
Equation (19) hence also reduces to (30) along the $\cal{C}^{+}$and
$\cal{C}^{-}$curves.
Finally, we consider the last differential form (18), expressing the
boundedness of $(\partial\psi_{2}/\partial y)$. Whether it is proportional to
(30) is not straightforward in particular due to the complex expression of
$C^{2}_{2y}$ and the ratio $C^{3}_{2y}/C^{2}_{2y}$ that is not obviously equal
to $t$. Nevertheless, we have proven that, along the characteristic curves the
differential form expressing the boundedness of $(\partial\psi_{2}/\partial
x)$ and $(\partial\psi_{2}/\partial y)$ (equations (15) and (18)) are
proportional by a minus $t$ factor. So we may use this property to derive a
simpler expression of $C^{2}_{2y}$ or prove the proportionality of (15) with
(30) along the $\cal{C}^{+}$and $\cal{C}^{-}$curves. Whatever the approach the
condition for proportionality reads

$t~{}(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)(u^{2}+v^{2})=(\gamma_{1}(ut+v)(E_{c}+H)+uv^{2}t-\gamma_{1}u^{2}v-\gamma
v^{3})(u+vt),$

that is found to be equivalent to (29) wherever $uv\neq 0$.
In a final formal calculation verification, it was checked that, on the
$\cal{C}^{+}$and $\cal{C}^{-}$curves, (17), is proportional to (30) (and we
already know that equations (18), (19), (20) are proportional to (17) along
these curves).

###  2.5 Main results and extension to 3D

We have found two differential equations, (24) and (25), valid along the
streamtraces $\cal{S}$for the adjoint system:

$E_{c}~{}\frac{d\psi_{1}}{ds}+H(u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds})+~{}H^{2}\frac{d\psi_{4}}{ds}=0\qquad\quad\frac{d\psi_{1}}{ds}+u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+E_{c}\frac{d\psi_{4}}{ds}=0.$

They are the counterpart of the constant total enthalpy and constant entropy
properties for the flow. They may be combined (for example to derive again
$\overline{U}.\nabla\psi_{1}-H\overline{U}.\nabla\psi_{4}=0$) or their
coefficients may be expressed differently using the well-known equations
satisfied by steady inviscid flows along a streamtrace. Let us finally note
that their straightforward 3D extensions (with natural notations),

$\displaystyle
E_{c}~{}\frac{d\psi_{1}}{ds}+H(u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+w~{}\frac{d\psi_{4}}{ds})+H^{2}~{}\frac{d\psi_{5}}{ds}=0$
$\displaystyle\frac{d\psi_{1}}{ds}+u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+w~{}\frac{d\psi_{4}}{ds}+E_{c}~{}\frac{d\psi_{5}}{ds}=0,$

are valid. This is the subject of Appendix C. The demonstration used in 3D
also provides another method for deriving the 2D equations.
We have found one differential equation for the $\cal{C}^{+}$and one
differential equation for the $\cal{C}^{-}$, equation (31) with relevant value
of $t^{\pm}$ for each curve (27)

$\displaystyle(u+vt^{\pm})~{}\frac{d\psi_{1}}{ds}+(u^{2}+v^{2})(\frac{d\psi_{2}}{ds}+t^{\pm}\frac{d\psi_{3}}{ds})+H(u+vt^{\pm})~{}\frac{d\psi_{4}}{ds}=0.$
$\displaystyle
t^{+}=\frac{uv+c\sqrt{u^{2}+v^{2}-c^{2}}}{u^{2}-c^{2}}\quad\textrm{for a
$\cal{C}^{+}$}\qquad\qquad
t^{-}=\frac{uv-c\sqrt{u^{2}+v^{2}-c^{2}}}{u^{2}-c^{2}}\quad\textrm{for a
$\cal{C}^{-}$}$

They are the counterpart of the differential forms satisfied by primitive flow
variables along the $\cal{C}^{+}$and $\cal{C}^{-}$. In the simpler cases,
where the classical angular relations (5) are valid, these relations may be
used to express differently the coefficients. We do not expect these relations
involving a 2D slope, $t$, to admit an extension in 3D.

##  3  Assessment of the adjoint ODEs

###  3.1 Consistency with known flow perturbation mechanisms

The adjoint vector is known to express the influence of a flow perturbation on
the associated QoI. Although discrete and continuous adjoint methods are
nowdays common tools for shape optimization, flow control and receptivity-
sensitivity-stability analysis, adjoint vectors are not easily interpreted.
The reason is that a single adjoint component, at a given location, is the
rate of the change in the QoI to the amplitude of a local perturbation in the
corresponding flow equation only (eg for the first component, mass injection
without perturbation of the momentum and energy equations). These individual
equation perturbations, of course, do not correspond to any realistic
possibility. A second complementary point of view, already mentioned in the
introduction, consists in calculating the dot product of the adjoint
components with the vector of a realistic perturbation and discuss the map of
the actuation influence [4, 9].
Plots of individual adjoint components for Euler flows appear in a 2002
publication by Venditti and Darmofal [19] (fig. 19 and 26). The discussion of
a x-momentum $CL$-adjoint plot in [19] mentions a singularity in the adjoint
along the stagnation streamline and a weak discontinuity upstream of the
primal shock on the upper surface. (With a finer mesh, the latter would have
been identified as a $\cal{C}^{-}$impinging the upperside shock-foot).
Although not discussed by the authors of [19], a corresponding plot for a
supersonic flow about two airfoils, strongly suggest backward information
propagation along the $\cal{S}$, $\cal{C}^{-}$and $\cal{C}^{+}$curves from the
support of the QoI. Concerning transonic airfoil flows, reference [20]
describes the mechanism by which locally perturbing one component of the flow
along the $\cal{C}^{-}$(resp. $\cal{C}^{+}$) impinging the upperside (resp.
lowerside) shock-foot results in a strong change in the lift and drag: the
flow perturbation propagates along the $\cal{C}^{-}$(resp. $\cal{C}^{+}$) and
results in a displacement of the shock.
The evaluation of the influence of physical source terms on the lift or drag
of a profile goes back to Giles and Pierce [4] who introduced the four
physical local source terms recalled in section I. Consistently with equations
(24), (25) and (31) the results obtained with this approach also support
adjoint information propagation along the $\cal{S}$, $\cal{C}^{-}$and
$\cal{C}^{+}$, from the support of the QoI (backwards w.r.t. the direction of
flow perturbations). Regarding the specific goal of acting at a shockfoot and
the four aforementioned source terms [4], the authors of [9] demonstrate that
the first two source terms (mass source at stagnation conditions, and normal
force) may displace the shock and strongly alter near-field forces if located
along the $\cal{C}^{+}$/ $\cal{C}^{-}$of the shock-foot whereas the fourth
source (change in stagnation pressure at fixed static pressure and total
enthalpy) may displace the shock-foot if located along the stagnation
streamline or along the wall upwind the shock.
The demonstrated ODEs along the $\cal{S}$, $\cal{C}^{-}$and
$\cal{C}^{+}$curves are hence consistent with known lines of specific
influence on drag or lift of classical steady Eulerian flows. These lines also
appear in the search of optimal forcings in control studies [21, 22] but we do
not extend on this aspect due to the different base equations.

###  3.2 Consistency with the analytical adjoint field of 2D supersonic
constant flow areas and the equations for the adjoint gradient at shocks

In reference [20] Todarello et al. derived the mathematical expression of the
2D Eulerian adjoint vector in a supersonic zone with constant flow (typically
upwind the detached shockwave created by an airfoil). The angle of attack
being $\alpha$ and Mach number $M$, this formula reads

$\displaystyle\psi(x,y)$ $\displaystyle=$
$\displaystyle\varphi_{\alpha}(x~{}\sin(\alpha)-y~{}\cos(\alpha))\lambda^{\alpha}_{0}+\varphi_{\alpha+\beta}(x~{}\sin(\alpha+\beta)-y~{}\cos(\alpha+\beta))\lambda^{\alpha+\beta}_{0}$
(32) $\displaystyle+$
$\displaystyle\varphi_{\alpha-\beta}(x~{}\sin(\alpha-\beta)-y~{}\cos(\alpha-\beta))\lambda^{\alpha-\beta}_{0},$

where $\beta=\textrm{sin}^{-1}(1/M)$, $\varphi_{\alpha}$,
$\varphi_{\alpha-\beta}$, $\varphi_{\alpha-\beta}$ are three scalar functions,
the $\lambda_{0}^{\mu}$ are left eigenvectors[9] of $A\sin(\mu)-B\cos(\mu)$

$\lambda_{0}^{\alpha-\beta}=\begin{pmatrix}\frac{c}{\rho}\left(1+\frac{\gamma_{1}}{2}M^{2}\right)\\\
\frac{1}{\rho}\left(\sin(\alpha-\beta)-\gamma_{1}\displaystyle\frac{u}{c}\right)\\\
\frac{1}{\rho}\left(-\cos(\alpha-\beta)-\gamma_{1}\displaystyle\frac{v}{c}\right)\\\
\frac{\gamma_{1}}{\rho
c}\end{pmatrix},\quad\lambda_{0}^{\alpha+\beta}=\begin{pmatrix}\frac{c}{\rho}\left(1+\frac{\gamma_{1}}{2}M^{2}\right)\\\
-\frac{1}{\rho}\left(\sin(\alpha+\beta)+\gamma_{1}\displaystyle\frac{u}{c}\right)\\\
-\frac{1}{\rho}\left(-\cos(\alpha+\beta)+\gamma_{1}\displaystyle\frac{v}{c}\right)\\\
\frac{\gamma_{1}}{\rho
c}\end{pmatrix}\quad\lambda_{0}^{\alpha}=\begin{pmatrix}-1-\frac{\gamma_{1}}{2}M^{2}\\\
\frac{\gamma_{1}u}{c^{2}}+\frac{2~{}\cos(\alpha)}{\cos(\alpha)u+\sin(\alpha)v}\\\
\frac{\gamma_{1}v}{c^{2}}+\frac{2~{}\sin(\alpha)}{\cos(\alpha)u+\sin(\alpha)v}\\\
-\frac{\gamma_{1}}{c^{2}}.\end{pmatrix}$ (33)

which formulas [23] have been simplified here using the null eigenvalues
relations valid in this specific context:
$u~{}\sin(\alpha)-v~{}\cos(\alpha)=0$, $\quad
u~{}\sin(\alpha-\beta)-v~{}\cos(\alpha-\beta)+c=0$, $\quad
u~{}\sin(\alpha+\beta)-v~{}\cos(\alpha+\beta)-c=0$.  Equation (32) is the
mathematical formula for the three stripes field depicted in figure 2. Each
$\varphi$ scalar function expresses the combined variation in the amplitude of
the adjoint components normally to the stripe direction. Each stripe is
crossed by the characteristic lines oriented in the direction of the other two
and we question whether equation (24), (25), (31) provide new information on
the $\varphi$ functions.
Considering the stripe oriented in the $\alpha-\beta$ direction, we first note
that equation (31) with the $t^{-}$ value, is automatically satisfied in its
geometrical domain since
$\varphi_{\alpha-\beta}(x~{}\sin(\alpha-\beta)-y~{}\cos(\alpha-\beta))\lambda^{\alpha-\beta}_{0}$
induces no variation of the $\psi$ components in the $\cal{C}^{-}$direction.
Regarding the conditions for satisfying (24) and (25) along an $\cal{S}$curve,
and satisfying (31) along a $\cal{C}^{+}$curve where they cross the
$\alpha-\beta$ stripe (as in fig. 2 right), we introduce the three functions

$\displaystyle\Gamma^{1}_{S}(s)$ $\displaystyle=$ $\displaystyle
E_{c}\psi_{1}+Hu\psi_{2}+Hv\psi_{3}+H^{2}\psi_{4}$
$\displaystyle\Gamma^{2}_{S}(s)$ $\displaystyle=$
$\displaystyle\psi_{1}+u\psi_{2}+v\psi_{3}+E_{c}\psi_{4}$
$\displaystyle\Gamma_{C+}(s)$ $\displaystyle=$
$\displaystyle(u+vt^{+})\psi_{1}+(u^{2}+v^{2})(\psi_{2}+t^{+}\psi_{3})+H(u+vt^{+})\psi_{4},$

with $s$ the curvilinear abscissa along the curve mentioned in the index. It
is proven in appendix D that

$\frac{d\Gamma^{1}_{S}}{ds}=0\qquad\frac{d\Gamma^{2}_{S}}{ds}=0\qquad\frac{d\Gamma_{C+}}{ds}=0$

without condition on the $\varphi$ functions. As the flow is constant, this is
equivalent to the satisfaction of the ODEs along the $\cal{S}$ and the
$\cal{C}^{+}$curves. The other possible crossing of the stripes and the curves
also do not provide conditions on the $\varphi$ functions. We hence have not
gained information on the analytical adjoint field of 2D supersonic constant
flow zones but this automatic consistency establishes a link, via a series the
orthogonality properties, between the coefficients of the differential
equation (24), (25), (31) and the relevant left eigenvectors of the Euler flux
Jacobian (33).

Regarding shock-waves, it is know that, for classical QoIs like pressure
integrals at the wall, the adjoint vector is continuous at shocks although the
flow is not, whereas the adjoint gradient may be discontinuous. In most common
cases, although a shockwave only makes the normal component of the velocity
subsonic, the $\cal{C}^{+}$and $\cal{C}^{-}$curves end at the shock and the
most interesting discussion regards the consequences of (24), (25). As they
are valid both sides of the shock, the jump operator may be applied to them
across the discontinuity $\Sigma$:

$\llbracket
E_{c}~{}\frac{d\psi_{1}}{ds}+H(u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds})+~{}H^{2}\frac{d\psi_{4}}{ds}\rrbracket=0\qquad\llbracket\frac{d\psi_{1}}{ds}+u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+E_{c}\frac{d\psi_{4}}{ds}\rrbracket=0,$

where all terms but $H$ may be discontinuous. Unfortunately, it does not seem
possible to reduce these equations to the simpler ones for the adjoint
gradient discontinuity [8, 9] where the derivatives in the two directions of
the local frame of reference attached to the shock appear independently.
Conversely, we note that in the case of a normal shock, the equations derived
in [8] prove the previous two jump equations.

Figure 2: Left: Sketch of lift/drag adjoint field for a supersonic flow about
an airfoil. Individual stripes (blue) and superimposition of the three stripes
(grey) in the constant flow zone ; nonzero adjoint zone downwind the shockwave
(green). Right: $\cal{C}^{-}$stripe crossed by a $\cal{S}$and a
$\cal{C}^{+}$curve

###  3.3 Numerical assessment method

The numerical assessment method consists in computing flow and adjoint fields
over a very fine mesh, and calculating the following integrals:

$\displaystyle K_{O}S^{1}$ $\displaystyle=$ $\displaystyle\int_{\cal
S}\left(E_{c}~{}\frac{d\psi_{1}}{ds}+H(u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds})+~{}H^{2}\frac{d\psi_{4}}{ds}\right)ds$
(34) $\displaystyle K_{O}S^{2}$ $\displaystyle=$ $\displaystyle\int_{\cal
S}\left(\frac{d\psi_{1}}{ds}+u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+~{}E_{c}\frac{d\psi_{4}}{ds}\right)ds$
(35) $\displaystyle K_{O}C^{+}$ $\displaystyle=$ $\displaystyle\int_{{\cal
C}^{+}}\left((u+vt^{+})~{}\frac{d\psi_{1}}{ds}+(u^{2}+v^{2})(\frac{d\psi_{2}}{ds}+t^{+}\frac{d\psi_{3}}{ds})+H(u+vt^{+})~{}\frac{d\psi_{4}}{ds}\right)ds$
(36) $\displaystyle K_{O}C^{-}$ $\displaystyle=$ $\displaystyle\int_{{\cal
C}^{-}}\left((u+vt^{-})~{}\frac{d\psi_{1}}{ds}+(u^{2}+v^{2})(\frac{d\psi_{2}}{ds}+t^{-}\frac{d\psi_{3}}{ds})+H(u+vt^{-})~{}\frac{d\psi_{4}}{ds}\right)ds.$
(37)

Here the intermediate subscript $O$ stands for the output functional of
interest ; it is subsequently replaced by $L$ (for the lift, $CLp$) and $D$
(for drag, $CDp$, consistently with the adjoint vector placed on the right-
hand side.

The integration is performed in the forward sense for the adjoint, that is,
backwards w.r.t. the direction of the flow information propagation. The
integration domain for the above line integrals extends to the interior of the
disk of radius $3c$ centred at $(0.5c,0)$, chosen for plotting readability,
while the flow computational domain itself extends to 150$c$. It may be
shorter, in particular in the transonic case where the $\cal{C}^{+}$and
$\cal{C}^{-}$curves are limited to the supersonic bubble(s). The four
quantities $K_{O}S^{1}$,…, $K_{O}C^{-}$ are expected to be close to zero and,
to avoid any error in scale, also calculated and plotted are the corresponding
subparts, that is, for $K_{O}S^{1}$ for example,

$K_{O}S^{1}_{1}=\int_{\cal S}\left(E_{c}~{}\frac{d\psi_{1}}{ds}\right)ds\qquad
K_{O}S^{1}_{2}=\int_{\cal S}\left(Hu~{}\frac{d\psi_{2}}{ds}\right)ds$
$K_{O}S^{1}_{3}=\int_{\cal S}\left(Hv~{}\frac{d\psi_{3}}{ds}\right)ds\qquad
K_{O}S^{1}_{4}=\int_{\cal S}\left(H^{2}\frac{d\psi_{4}}{ds}\right)ds.$

The sum of the four terms is expected to be much smaller than each one of them
individually. All the integrals are calculated backwards, along a finely
discretized characteristic curve, simply by the trapezoidal rule.
The discrete flows and adjoints were available from former computations [9] in
which the Jameson-Schmidt-Turkel scheme [24] was applied, and using the
discrete adjoint module of the $elsA$ code [25, 26]. Of course when trying to
assess properties of exact adjoint fields from numerical discrete solutions,
it is desirable to work either with continuous adjoints or with dual
consistent discrete adjoints[27, 28, 29]. Precisely in [9], it was
demonstrated for structured meshes how to slightly modify the scheme’s
Jacobian (in the derivative of the dissipation flux, for the next to wall
faces) to get a dual consistent linearization. This slight modification of the
exact scheme Jacobian is retained here to work with adjoint fields that are
consistent with the continuous equations discussed in §2. Note also that these
adjoint fields have also been satisfactorily verified by a posteriori
discretization of the continuous adjoint equation[9].
Only the solutions calculated over the finest mesh defined in reference [30]
(structured 4097$\times$4097 mesh) are used here. The iso-Mach number lines,
iso-first component of $CLp$ adjoint and the extracted curves may be seen for
all cases in figure 3.

Figure 3: Simulations over the 4097$\times$4097 mesh of [30]. iso-Mach number
lines (upper three plots) and iso-$\psi^{1}_{CLp}$ lines (lower three plots)
and extracted curves. Left, subsonic case: streamline. Middle, transonic case:
streamline, $\cal{C}^{+}$, $\cal{C}^{-}$(black) and sonic line (violet). Right
supersonic case: streamline, $\cal{C}^{+}$(upper side), $\cal{C}^{-}$(lower
side)

###  3.4 Numerical assessment of the ODEs for a supersonic flow about the
NACA0012 airfoil

The retained flow conditions are $M_{\infty}=1.5~{},~{}\alpha=1^{o}$. We first
assess the streamtraces equations (24) (25). The $K_{D}S^{1}$, $K_{D}S^{2}$,
$K_{L}S^{1}$, $K_{L}S^{2}$ integrals and their subparts are calculated along
the trajectory passing through $(c,0.1c)$. The integration indeed leads to
very small values of $K_{D}S^{1}$, $K_{D}S^{2}$, $K_{L}S^{1}$, $K_{L}S^{2}$
along the curve w.r.t. their subparts. It is well-known for this of kind of
flow that the exact lift- and drag- adjoint is equal to zero downstream the
backwards flow-characteristics emanating from the trailing edge (since no
perturbation downstream those two lines can affect the pressure on the
aerofoil and, consequently, the lift or the drag – see for example fig. 6 and
A21 in [9]). This property is well satisfied by discrete adjoint fields and,
as the integration is performed backwards along the streamtrace, null values
of $K_{D}S^{1}$, $K_{D}S^{2}$, $K_{L}S^{1}$, $K_{L}S^{2}$ and all their
subparts are observed above a specific $x$ corresponding to the intersection
of the streamtrace with this trailing-edge $\cal{C}^{-}$. The integration of
(24) and (25) reveals (i) the discontinuity of the integration variables (the
adjoint components) at $x\simeq 0.85$ that appears as a discontinuity in the
subpart curves ; (ii) a discontinuity of the integration coefficients, when
crossing the detached shock wave at $x\simeq-0.08$, that results in strong
gradients in the subparts curves. None of those discontinuities alters the
almost null values of $K_{D}S^{1}$ $K_{D}S^{2}$, $K_{D}L^{1}$ and
$K_{L}S^{2}$. The strong adjoint gradients in the subsonic bubble
(approximately $x\in~{}[-0.08,0.05]$) also clearly translate the $K_{S}$
curves. Finally, note that for $x$ lower than $\simeq-0.05$ the flow is
constant (so the streamtrace is a straight line) and for $x$ lower than
$\simeq-1$ the backward streamtrace enters a zone of constant adjoint along
its direction (see fig. 2 left and equation (32)). This theoretical property
is well translated in constant sections of the curves. The results are
equivalently accurate for lift and drag. They are presented for the drag in
figure 4.

Figure 4: $M_{\infty}=1.50$, $\alpha=1^{o}$, (4097$\times$4097 mesh [30])
Numerical assessment of equation (24) (left) and (25) (right) for the lift.
Method of verification: the black curve should ideally coincide with the $x$
axis

A $\cal{C}^{+}$and a $\cal{C}^{-}$curve are then extracted using equation (5).
The selected $\cal{C}^{+}$is initiated at $x\simeq 0.3$ upper side and the
retained $\cal{C}^{-}$starts at the same abscissa but on the lower side. This
choice was guided by the extraction method and the observation that $k^{+}$
(resp. $k^{-}$) is almost constant on the lower (resp. upper) side. The
$K_{D}C^{+}$, $K_{L}C^{+}$, $K_{D}C^{-}$,$K_{L}C^{-}$ terms and their subparts
have been computed. The results appear to be very satisfactory. Also observed
is the equality of $K_{D}C^{+}_{1}$ and $K_{D}C^{+}_{4}$, $K_{L}C^{+}_{1}$ and
$K_{L}C^{+}_{4}$ and along the $\cal{C}^{+}$and correspondingly along the
$\cal{C}^{-}$curves. This is due to the fact that $\psi_{1}=H\psi_{4}$ (for
Euler flows and for pressure-based integrals along the wall that is well
satisfied at the discrete level – see for example [9] fig. A21, A22, A23) and
to the expression of the $d\psi_{1}$ and $d\psi_{4}$ terms in (30). Figure 5
presents the verification plots for the two functions along the selected
$\cal{C}^{+}$.

Figure 5: $M_{\infty}=1.50$, $\alpha=1^{o}$, (4097$\times$4097 mesh [30])
Numerical assessment of equation (30) for the lift (left) and for the drag
(right). Method of verification: the black curve should ideally coincide with
the $x$ axis

###  3.5 Numerical assessment of the ODEs for a transonic flow about the
NACA0012 airfoil

The flow conditions are $M_{\infty}=0.85~{},~{}\alpha=2^{o}$. Careful
verification of the streamstraces equations (24) (25) has been performed for
the streamtrace passing through $(c,0.1c)$ and very satisfactory results have
been found. As in the previous section, the intersection of the $\cal{S}$curve
with the shockwave results in the subsparts curves in very strong gradients
but does not affect the almost null value of $K_{D}S^{1}$, $K_{D}S^{2}$,
$K_{L}S^{1}$, $K_{L}S^{2}$. As similar results have been presented in the
previous and the next subsections, we focus here on the $\cal{C}^{+}$and
$\cal{C}^{-}$curves. A $\cal{C}^{+}$and a $\cal{C}^{-}$curves have been
extracted taking care to select the longest possible curves (and to avoid, for
the $\cal{C}^{-}$the curve passing by the shock-foot where numerical
divergence of the adjoint may be observed). The selected $\cal{C}^{+}$(resp.
$\cal{C}^{-}$) is passing approximately through the point $(0.197,0.057)$
(resp. $(0.954,0141)$). The verification of the consistency of the numerical
solutions w.r.t. (30) is satisfactory although the largest observed errors
appear in this case, for the $\cal{C}^{-}$curve, for the lift, in the vicinity
of the inlet of the supersonic bubble – see figure 6. This largest observed
error is about 2% of the largest absolute value of the four subparts. The
integrals along the $\cal{C}^{-}$are regular, whereas those along the
$\cal{C}^{+}$exhibit a sharp peak close to $x\simeq 0.52$, at the intersection
with the $\cal{C}^{-}$passing by the shockfoot. We refer to §3.1 for the
reason of the corresponding strong adjoint gradient.

Figure 6: $M_{\infty}=0.85$, $\alpha=2^{o}$, (4097$\times$4097 mesh [30])
Numerical assessment of equation (30) for the lift (left) and for the drag
(right), for the selected $\cal{C}^{-}$(up) and $\cal{C}^{-}$(down). Method of
verification: the black curve should ideally coincide with the $x$ axis

###  3.6 Numerical assessment of the streamtrace ODEs for a subsonic flow
about the NACA0012 airfoil

We expect relations (24) and (25) to be valid along the trajectories of a
subsonic flow. The retained flow conditions have been
$M_{\infty}=0.4~{},~{}\alpha=5^{o}$. The $K_{D}S^{1}$, $K_{D}S^{2}$,
$K_{L}S^{1}$, $K_{L}S^{2}$ integrals and their subparts are calculated along
the trajectory passing through $(c,0.1c)$. The integration indeed leads to
very small values of $K_{D}S^{1}$, $K_{D}S^{2}$, $K_{L}S^{1}$, $K_{L}S^{2}$
along the curve w.r.t. their subparts. The results are equivalently good for
lift and drag and are presented in figure 7 for the lift. The lower left plot
of figure 3 presents the typical aspect of subsonic lift or drag adjoints
indicating that an actuation able to significantly alter these QoIs is to be
applied in the immediate vicinity of the wall. This property translates in
weakly varying curves far from the profile in our verification plots for this
test case.

Figure 7: $M_{\infty}=.40$, $\alpha=5^{o}$, (4097$\times$4097 mesh [30])
Numerical assessment of equation (24) (left) and (25) (right) for the lift.
Method of verification: black curve should ideally coincide with the $x$ axis

##  4 Conclusion

Ordinary Differential Equations have been derived for the adjoint Euler
equations using in a first step the method of characteristics in 2D. The
differential equations satisfied along the streamtraces in 2D have then been
extended to 3D and the combination of equations method used for the derivation
in this case also provides a simpler proof of the corresponding 2D equations.
All these ODEs are non linear differential equations that cannot be integrated
for non-constant flows.
The adjoint vector expresses the sensitivity of its corresponding scalar QoI
to a local perturbation in the flow equations. Its variations in the fluid
domain are often difficult to analyse as it precisely avoids the calculation
of the flow perturbation that causes the change in the QoI. Nevertheless lift
and drag adjoint fields have been examined since a long time, and the
presented equations for 2D problems clarify their highest values and the
strong sensitivity of their associated QoIs to perturbations located along
specific lines [19, 20, 8, 9].
These findings have been illustrated with flows, lift-adjoints, and drag-
adjoints over the classical NACA0012 airfoil using a very fine mesh and a
dual-consistent adjoint method. The conducted tests lead to very satisfactory
results (although minor deviations in our transonic case close to the inlet of
the upper side supersonic bubble). The demonstrated equations (24),(25), and
(31) hence also provide a verification tool for discrete adjoint fields.

SUPPLEMENTARY MATERIAL
Supplementary material consists of:
– five python scripts allowing to check the algebraic expressions of the
$C^{l}_{mx}$ and $C^{l}_{my}$ coefficients w.r.t. their definition as
determinants and one python script allowing the numerical verification of the
results of IIIB ;
– six outputs of Maple scripts checking the rank of the sets of differential
forms satisfied along the streamtraces, $\cal{C}^{+}$, and $\cal{C}^{-}$.

ACKNOWLEDGMENTS
The authors express their warm gratitude to J.C. Vassberg and A. Jameson for
allowing the co-workers of D. Destarac to use their hierarchy of O-grids
around the NACA0012 airfoil, as well as E. Hubert and B. Mourrain (Inria
Aromath Project Team) for their helpful guidance for the analytical
developments using the computer algebra software Maple. The authors also thank
Florent Renac, Fulvio Sartor and Martin Duguey for many fruitful discussions.

This research did not receive any specific grant from funding agencies in the
public, commercial or not-for-profit sector.

DATA AVAILABILITY
The data consists of the three flow-fields about the NACA0012 airfoil and the
corresponding lift-adjoint and drag-adjoint. They are available as Tecplot
binary or formatted files from the corresponding author.

AUTHOR CONTRIBUTIONS
Jacques Peter: Investigation (lead); Methodology (lead); Validation (equal);
Writing original draft (lead). Jean-Antoine Désidéri: Investigation
(supporting); Methodology (supporting); Validation (equal); Writing original
draft (supporting); Writing review and editing.

##  Appendix A Calculation and expressions of the $C_{ij}^{k}$ coefficients

Two or three of the vectors $(-B_{1}+tA_{1})$, $(-B_{2}+tA_{2})$,
$(-B_{3}+tA_{3})$ and $(-B_{4}+tA_{4})$ appear in the formulas of the
$C_{ij}^{k}$ coefficients expressed as the determinant of a $4\times 4$
matrix. They may be precalculated as

$-B_{1}+tA_{1}=\begin{bmatrix}0\\\ t\\\ -1\\\ 0\end{bmatrix}\qquad-
B_{2}+tA_{2}=\begin{bmatrix}t~{}\gamma_{1}~{}E_{c}-u\kappa\\\
-v+t(3-\gamma)u\\\ -u-t\gamma_{1}v\\\ t\gamma_{1}\end{bmatrix}$
$-B_{3}+tA_{3}=\begin{bmatrix}-\gamma_{1}~{}E_{c}-v\kappa\\\ \gamma_{1}u+tv\\\
-(3-\gamma)v+tu\\\ -\gamma_{1}\end{bmatrix}\qquad-
B_{4}+tA_{4}=\begin{bmatrix}(\gamma_{1}~{}E_{c}-H)\kappa\\\
tH-\gamma_{1}u\kappa\\\ -H-\gamma_{1}v\kappa\\\ \gamma\kappa\end{bmatrix}$

The expressions of the $C^{l}_{1x}$, $C^{l}_{2x}$, $C^{l}_{3x}$, $C^{l}_{1y}$,
$C^{l}_{2y}$, $C^{l}_{3y}$, $C^{l}_{4y}$ are gathered below (the $C^{l}_{4x}$
being given in §2.2). We recall that, with our notations, the null
differential form along the $\cal{S}$, $\cal{C}^{-}$and $\cal{C}^{+}$derived
from the existence of, e.g. $(\partial\psi_{i}/\partial y)$, reads
$C_{iy}^{1}d\psi_{1}-C_{iy}^{2}d\psi_{2}+C_{iy}^{3}d\psi_{3}-C_{iy}^{4}d\psi_{4}=0$

$\bullet$ The coefficients $C^{l}_{1x}$ are expressed below

$\displaystyle C^{1}_{1x}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\kappa~{}((2tu+(t^{2}-1)v)(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)-(u+vt)(\gamma_{1}tH+\gamma_{1}u\kappa+\gamma v\kappa t))$
$\displaystyle C^{2}_{1x}$ $\displaystyle=$ $\displaystyle
dx^{2}dy~{}\gamma_{1}~{}\kappa~{}(u^{2}+v^{2})~{}H$ $\displaystyle C^{3}_{1x}$
$\displaystyle=$ $\displaystyle-
dx^{2}dy~{}\gamma_{1}~{}\kappa~{}t~{}(u^{2}+v^{2})~{}H$ $\displaystyle
C^{4}_{1x}$ $\displaystyle=$ $\displaystyle
dx^{2}dy~{}\gamma_{1}~{}\kappa~{}(u+vt)~{}H^{2}$

$\bullet$ The coefficients $C^{l}_{2x}$ are expressed below

$\displaystyle C^{1}_{2x}$ $\displaystyle=$ $\displaystyle-
dx^{2}dy~{}\kappa~{}(\gamma_{1}H+\gamma_{1}E_{c}+\gamma v\kappa)$
$\displaystyle C^{2}_{2x}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\kappa~{}(\gamma_{1}u^{2}v-uv^{2}t+\gamma
v^{3}-\gamma_{1}(ut+v)(E_{c}+H))$ $\displaystyle C^{3}_{2x}$ $\displaystyle=$
$\displaystyle dx^{2}dy~{}\kappa~{}(\gamma_{1}u^{2}v-uv^{2}t+\gamma
v^{3}-\gamma_{1}(ut+v)(E_{c}+H))$ $\displaystyle C^{4}_{2x}$ $\displaystyle=$
$\displaystyle dx^{2}dy~{}H~{}\kappa~{}(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)$

$\bullet$ The coefficients $C^{l}_{3x}$ are expressed below

$\displaystyle C^{1}_{3x}$ $\displaystyle=$ $\displaystyle
dx^{2}dy~{}\kappa~{}(t\gamma_{1}H+t\gamma_{1}E_{c}-\gamma u\kappa)$
$\displaystyle C^{2}_{3x}$ $\displaystyle=$ $\displaystyle-
dx^{2}dy~{}\kappa~{}((\gamma+1)u^{2}v-t\gamma_{1}uv^{2}+(v+ut)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2}))$ $\displaystyle C^{3}_{3x}$ $\displaystyle=$ $\displaystyle
dx^{3}~{}\kappa~{}((tv^{2}-u\kappa)(-\gamma_{1}u+(\gamma+1)tv)-(ut-(1+2t^{2})v)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
v^{2}))$ $\displaystyle C^{4}_{3x}$ $\displaystyle=$ $\displaystyle-
dx^{2}dy~{}H~{}\kappa~{}(t\gamma_{1}H+t\gamma_{1}E_{c}-\gamma u\kappa)$

$\bullet$ The coefficients $C^{l}_{1y}$ are expressed below

$\displaystyle C^{1}_{1y}$ $\displaystyle=$ $\displaystyle
dx^{3}~{}\kappa~{}\left((\gamma_{1}+\gamma
t^{2})u^{3}-2u^{2}vt-\gamma_{1}t\kappa
H+(\gamma+\gamma_{1}t^{2})uv^{2}-\gamma_{1}u~{}(1+t^{2})~{}E_{c}\right)$
$\displaystyle C^{2}_{1y}$ $\displaystyle=$ $\displaystyle-
dx^{3}\gamma_{1}~{}\kappa~{}(u^{2}+v^{2})~{}H$ $\displaystyle C^{3}_{1y}$
$\displaystyle=$ $\displaystyle
dx^{3}~{}\gamma_{1}~{}t~{}\kappa~{}(u^{2}+v^{2})~{}H$ $\displaystyle
C^{4}_{1y}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\gamma_{1}~{}\kappa~{}(u+vt)~{}H^{2}$

$\bullet$ The coefficients $C^{l}_{2y}$ are expressed below

$\displaystyle C^{1}_{2y}$ $\displaystyle=$ $\displaystyle
dx^{3}~{}\kappa~{}(\gamma_{1}H+\gamma_{1}E_{c}+\gamma v\kappa)$ $\displaystyle
C^{2}_{2y}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\kappa~{}((v\kappa+u^{2})(vt-(\gamma_{1}+\gamma
t^{2})u)-(vt-(2+t^{2})u))(\gamma_{1}E_{c}+\gamma_{1}H+\gamma v\kappa))$
$\displaystyle C^{3}_{2y}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\kappa~{}(\gamma_{1}u^{2}v-uv^{2}t+\gamma
v^{3}-\gamma_{1}(ut+v)(E_{c}+H))$ $\displaystyle C^{4}_{2y}$ $\displaystyle=$
$\displaystyle-dx^{3}~{}H~{}\kappa~{}(\gamma_{1}H+\gamma_{1}E_{c}+\gamma
v\kappa)$

$\bullet$ The coefficients $C^{l}_{3y}$ are expressed below

$\displaystyle C_{3y}^{1}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\kappa~{}(\gamma_{1}tE_{c}+\gamma_{1}tH-\gamma u\kappa)$
$\displaystyle C_{3y}^{2}$ $\displaystyle=$ $\displaystyle
dx^{3}~{}\kappa~{}((\gamma+1)u^{2}v-t\gamma_{1}uv^{2}+(v+ut)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2}))$ $\displaystyle C_{3y}^{3}$ $\displaystyle=$ $\displaystyle-
dx^{3}t\kappa((\gamma+1)u^{2}v-t\gamma_{1}uv^{2}+(v+ut)(\gamma_{1}E_{c}+\gamma_{1}H-\gamma
u^{2}))$ $\displaystyle C_{3y}^{4}$ $\displaystyle=$ $\displaystyle
dx^{3}\kappa~{}H~{}\left(\gamma_{1}tE_{c}+\gamma_{1}tH-\gamma u\kappa\right)$

$\bullet$ The coefficients $C^{l}_{4y}$ are expressed below

$\displaystyle C_{4y}^{1}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\gamma_{1}~{}\kappa~{}(u+vt)$ $\displaystyle C_{4y}^{2}$
$\displaystyle=$ $\displaystyle dx^{3}~{}\gamma_{1}~{}\kappa~{}(u^{2}+v^{2})$
$\displaystyle C_{4y}^{3}$ $\displaystyle=$ $\displaystyle-
dx^{3}~{}\gamma_{1}~{}\kappa~{}t~{}(u^{2}+v^{2})$ $\displaystyle C_{4y}^{4}$
$\displaystyle=$ $\displaystyle dx^{3}\kappa\left((\gamma_{1}+\gamma
t^{2})u^{3}-2u^{2}vt-\gamma_{1}t\kappa
H+(\gamma+\gamma_{1}t^{2})uv^{2}-\gamma_{1}u~{}(1+t^{2})~{}E_{c}\right)$

##  Appendix B General properties of the $C^{l}_{mx}$ and $C^{l}_{my}$
coefficients

Equation (7) refers to the limit of small space steps and the search of
characteristic curves ; nevertheless, the expressions of the $C^{l}_{mx}$ and
$C^{l}_{my}$ coefficients may be considered for an arbitrary direction and an
arbitrary norm of vector $(dx,dy)$. Without any assumption linking $(dx,dy)$
and $(u,v)$, the relations between the coefficients of the same differential
forms are:

$\displaystyle C^{4}_{1x}$ $\displaystyle=$ $\displaystyle HC^{1}_{1x}\quad
C^{3}_{1x}=tC^{2}_{1x}~{}~{};~{}C^{4}_{2x}=-HC^{1}_{2x}\quad
C^{3}_{2x}=-tC^{2}_{2x}~{}~{};~{}C^{4}_{3x}=-HC^{1}_{3x}~{}~{};~{}C^{3}_{4x}=-tC^{2}_{4x}\qquad$
(38) $\displaystyle C^{4}_{1y}$ $\displaystyle=$ $\displaystyle
HC^{1}_{1y}\quad
C^{3}_{1y}=tC^{2}_{1y}~{}~{};~{}C^{4}_{2y}=-HC^{1}_{2y}~{}~{};~{}~{}C^{4}_{3y}=-HC^{1}_{3y}\quad
C^{3}_{3y}=-tC^{2}_{3y}~{}~{};~{}C^{3}_{4y}=-tC^{2}_{4y}\qquad$ (39)

Besides, twelve of the sixteen coefficients of the differential forms for the
$x$ and $y$ derivatives are proportional by a $(-t)$ factor:

$\displaystyle C^{2}_{1x}$ $\displaystyle=$
$\displaystyle-t~{}C^{2}_{1y}\qquad\quad C^{3}_{1x}=-t~{}C^{3}_{1y}\qquad\quad
C^{4}_{1x}=-t~{}C^{4}_{1y}\qquad\quad$ (40) $\displaystyle C^{1}_{2x}$
$\displaystyle=$ $\displaystyle-t~{}C^{1}_{2y}\qquad\quad
C^{3}_{2x}=-t~{}C^{3}_{2y}\qquad\quad C^{4}_{2x}=-t~{}C^{4}_{2y}\qquad\quad$
(41) $\displaystyle C^{1}_{3x}$ $\displaystyle=$
$\displaystyle-t~{}C^{1}_{3y}\qquad\quad C^{2}_{3x}=-t~{}C^{2}_{3y}\qquad\quad
C^{4}_{3x}=-t~{}C^{4}_{3y}\qquad\quad$ (42) $\displaystyle C^{1}_{4x}$
$\displaystyle=$ $\displaystyle-t~{}C^{1}_{4y}\qquad\quad
C^{2}_{4x}=-t~{}C^{2}_{4y}\qquad\quad C^{3}_{4x}=-t~{}C^{3}_{4y}\qquad\quad$
(43)

Finally, the $C^{1}_{1l}$ and $C^{4}_{4l}$ coefficients are equal

$C^{1}_{1x}=C^{4}_{4x}\qquad\qquad\qquad C^{1}_{1y}=C^{4}_{4y}$ (44)

##  Appendix C Streamstrace ODEs in dimension 3

From 2D equations (24) and (25), we can infer corresponding candidate
equations for the adjoint vector along the streamtraces in 3D:

$\displaystyle
E_{c}~{}\frac{d\psi_{1}}{ds}+H(u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+w~{}\frac{d\psi_{4}}{ds})+~{}H^{2}\frac{d\psi_{5}}{ds}=0$
(45)
$\displaystyle\frac{d\psi_{1}}{ds}+u~{}\frac{d\psi_{2}}{ds}+v~{}\frac{d\psi_{3}}{ds}+w~{}\frac{d\psi_{4}}{ds}+E_{c}\frac{d\psi_{5}}{ds}=0.$
(46)

Could these equations be possibly proven from the 3D Euler adjoint equations

$-A^{T}\frac{\partial\psi}{\partial x}-B^{T}\frac{\partial\psi}{\partial
y}-C^{T}\frac{\partial\psi}{\partial z}=0,$ (47)

where the 3D transposed Jacobian of Euler fluxes read

$A^{T}=\begin{bmatrix}0&(\gamma_{1}E_{c}-u^{2})&-uv&-uw&(\gamma_{1}E_{c}-H)u\\\
1&-\gamma_{1}u+2u&v&w&(H-\gamma_{1}u^{2})\\\
0&-\gamma_{1}v&u&0&-\gamma_{1}uv\\\ 0&-\gamma_{1}w&0&u&-\gamma_{1}uw\\\
0&\gamma_{1}&0&0&\gamma u\\\
\end{bmatrix}~{}~{}~{}~{}~{}~{}~{}~{}B^{T}=\begin{bmatrix}0&-uv&(\gamma_{1}E_{c}-v^{2})&-vw&(\gamma_{1}E_{c}-H)v\\\
0&v&-\gamma_{1}u&0&-\gamma_{1}uv\\\
1&u&-\gamma_{1}v+2v&w&(H-\gamma_{1}v^{2})\\\
0&0&-\gamma_{1}w&v&-\gamma_{1}vw\\\ 0&0&\gamma_{1}&0&\gamma v\\\
\end{bmatrix}$
$C^{T}=\begin{bmatrix}0&-uw&-vw&(\gamma_{1}E_{c}-w^{2})&(\gamma_{1}E_{c}-H)w\\\
0&w&0&-\gamma_{1}u&-\gamma_{1}uw\\\ 0&0&w&-\gamma_{1}v&-\gamma_{1}vw\\\
1&u&v&-\gamma_{1}w+2w&(H-\gamma_{1}w^{2})\\\ 0&0&0&\gamma_{1}&\gamma W\\\
\end{bmatrix}$

$s$ being the curvilinear abscissa along a trajectory in the direction
opposite to the flow displacement, the differentiation w.r.t. $s$ may be
expressed as

$\frac{d}{ds}=-\frac{u}{\|U\|}\frac{d}{dx}-\frac{v}{\|U\|}\frac{d}{dy}-\frac{w}{\|U\|}\frac{d}{dz}$

First considering (46), the equation with the simpler coefficients, this
equation is satisfied in the fluid domain if and only if

$\displaystyle(u\frac{d\psi_{1}}{dx}+v\frac{d\psi_{1}}{dy}+w\frac{d\psi_{1}}{dz})~{}+u(u\frac{d\psi_{2}}{dx}+v\frac{d\psi_{2}}{dy}+w\frac{d\psi_{2}}{dz})~{}$
$\displaystyle+v(u\frac{d\psi_{3}}{dx}+v\frac{d\psi_{3}}{dy}+w\frac{d\psi_{3}}{dz})~{}+w(u\frac{d\psi_{4}}{dx}+v\frac{d\psi_{4}}{dy}+w\frac{d\psi_{4}}{dz})~{}$
$\displaystyle+E_{c}(u\frac{d\psi_{5}}{dx}+v\frac{d\psi_{5}}{dy}+w\frac{d\psi_{5}}{dz})=0\
$ (48)

(where we have removed the norm of the velocity and the minus sign thanks to
the homogeneity of the equation). We now search if a combination of the lines
of (47) that would result in (48). For the required $\psi_{1}$ terms to
appear, the combination $(-uL_{2}-vL_{3}-wL_{4})$ ($L_{j}$ denoting the $j$-th
line of (47)) needs to be calculated:

$\displaystyle\left[u\right.$ $\displaystyle(-2\gamma_{1}E_{c}+2u^{2})$
$\displaystyle 2uv$ $\displaystyle 2uw$
$\displaystyle\left.(Hu-2\gamma_{1}uE_{c})\right]$
$\displaystyle\frac{\partial\psi}{\partial x}$ $\displaystyle+$
$\displaystyle\left[v\right.$ $\displaystyle 2uv$
$\displaystyle(-2\gamma_{1}E_{c}+2v^{2})$ $\displaystyle 2vw$
$\displaystyle\left.(Hv-2\gamma_{1}vE_{c})\right]$
$\displaystyle\frac{\partial\psi}{\partial y}$ $\displaystyle+$
$\displaystyle\left[w\right.$ $\displaystyle 2uw$ $\displaystyle 2vw$
$\displaystyle(-2\gamma_{1}E_{c}+2w^{2})$
$\displaystyle\left.(Hw-2\gamma_{1}wE_{c})\right]$
$\displaystyle\frac{\partial\psi}{\partial z}=0$

Subtracting the first line, calculating $(-L_{1}-uL_{2}-vL_{3}-wL_{4})$,
almost fixes the expected coefficients for the $\psi_{1}$ to $\psi_{4}$
derivatives:

$\displaystyle\left[u\right.$ $\displaystyle(-\gamma_{1}E_{c}+u^{2})$
$\displaystyle uv$ $\displaystyle uw$
$\displaystyle\left.-\gamma_{1}uE_{c}\right]$
$\displaystyle\frac{\partial\psi}{\partial x}$ $\displaystyle+$
$\displaystyle\left[v\right.$ $\displaystyle uv$
$\displaystyle(-\gamma_{1}E_{c}+v^{2})$ $\displaystyle vw$
$\displaystyle\left.-\gamma_{1}vE_{c}\right]$
$\displaystyle\frac{\partial\psi}{\partial y}$ $\displaystyle+$
$\displaystyle\left[w\right.$ $\displaystyle uw$ $\displaystyle vw$
$\displaystyle(-\gamma_{1}E_{c}+w^{2})$
$\displaystyle\left.-\gamma_{1}wE_{c}\right]$
$\displaystyle\frac{\partial\psi}{\partial z}=0$

Finally forming $(-L_{1}-uL_{2}-vL_{3}-wL_{4}-E_{c}L_{5})$ yields

$\displaystyle\left[u\right.$ $\displaystyle u^{2}$ $\displaystyle uv$
$\displaystyle uw$ $\displaystyle\left.uE_{c}\right]$
$\displaystyle\frac{\partial\psi}{\partial x}$ $\displaystyle+$
$\displaystyle\left[v\right.$ $\displaystyle uv$ $\displaystyle v^{2}$
$\displaystyle vw$ $\displaystyle\left.vE_{c}\right]$
$\displaystyle\frac{\partial\psi}{\partial y}$ $\displaystyle+$
$\displaystyle\left[w\right.$ $\displaystyle uw$ $\displaystyle vw$
$\displaystyle w^{2}$ $\displaystyle\left.wE_{c}\right]$
$\displaystyle\frac{\partial\psi}{\partial z}=0$

that is exactly equation (48).
In order to demonstrate (45), we first form the combination
$-((2E_{c}-H)L_{1}+uE_{c}L_{2}+vE_{c}L_{3}+wE_{c}L_{4})$ of the lines of the
3D adjoint Euler equations. This results in

$\displaystyle\left[uE_{c}\right.$ $\displaystyle(u^{2}H-\gamma_{1}HE_{c})$
$\displaystyle uvH$ $\displaystyle uw$ $\displaystyle\left.-\gamma
HuE_{c}+H^{2}u\right]$ $\displaystyle\frac{\partial\psi}{\partial x}$
$\displaystyle+$ $\displaystyle\left[vE_{c}\right.$ $\displaystyle uvH$
$\displaystyle(v^{2}H-\gamma_{1}HE_{c})$ $\displaystyle vw$
$\displaystyle\left.-\gamma HvE_{c}+H^{2}v\right]$
$\displaystyle\frac{\partial\psi}{\partial y}$ $\displaystyle+$
$\displaystyle\left[wE_{c}\right.$ $\displaystyle uwH$ $\displaystyle vwH$
$\displaystyle(w^{2}H-\gamma_{1}HE_{c})$ $\displaystyle\left.-\gamma
HwE_{c}+H^{2}w\right]$ $\displaystyle\frac{\partial\psi}{\partial z}=0$

Adding $-E_{c}~{}H~{}L_{5}$ yields

$\displaystyle\left[uE_{c}\right.$ $\displaystyle u^{2}H$ $\displaystyle uvH$
$\displaystyle uwH$ $\displaystyle\left.H^{2}u\right]$
$\displaystyle\frac{\partial\psi}{\partial x}$ $\displaystyle+$
$\displaystyle\left[vE_{c}\right.$ $\displaystyle uvH$ $\displaystyle v^{2}H$
$\displaystyle vwH$ $\displaystyle\left.H^{2}v\right]$
$\displaystyle\frac{\partial\psi}{\partial y}$ $\displaystyle+$
$\displaystyle\left[wE_{c}\right.$ $\displaystyle uwH$ $\displaystyle vwH$
$\displaystyle w^{2}H$ $\displaystyle\left.H^{2}w\right]$
$\displaystyle\frac{\partial\psi}{\partial z}=0$

that is equivalent to (45).

##  Appendix D Demonstration of the properties presented in §3.2

The differentiation of the functions $\Gamma^{1}_{S}$, $\Gamma^{2}_{S}$ and
$\Gamma_{C+}$ appearing in section §3.2 uses the expression of the adjoint
field where the $(\alpha-\beta)$-oriented stripe is non superimposed with the
two other:

$\psi(x,y)=\varphi_{\alpha-\beta}(x~{}\sin(\alpha-\beta)-y~{}\cos(\alpha-\beta))\lambda^{\alpha-\beta}_{0}$
$\displaystyle\frac{d\Gamma^{1}_{S}}{ds}$ $\displaystyle=$ $\displaystyle
E_{c}\frac{d\psi_{1}}{ds}+Hu\frac{d\psi_{2}}{ds}+Hv\frac{d\psi_{3}}{ds}+H^{2}\frac{d\psi_{4}}{ds}$
$\displaystyle\frac{d\Gamma^{1}_{S}}{ds}$ $\displaystyle=$
$\displaystyle\left(E_{c}(\frac{c}{\rho}(1+\frac{\gamma_{1}}{2}M^{2}))+Hu(\frac{1}{\rho}(\sin(\alpha-\beta)-\gamma_{1}\displaystyle\frac{u}{c}))+Hv(\frac{1}{\rho}(-\cos(\alpha-\beta)-\gamma_{1}\displaystyle\frac{v}{c}))+H^{2}(\frac{\gamma_{1}}{\rho
c})\right)\frac{d\varphi_{\alpha-\beta}}{d\xi}\frac{d\xi}{ds}$
$\displaystyle\frac{d\Gamma_{S}}{ds}$ $\displaystyle=$
$\displaystyle\left(\frac{\gamma_{1}E_{c}}{\rho
c}+\frac{c^{2}u}{\gamma_{1}\rho}\sin(\alpha-\beta)-\frac{c^{2}v}{\gamma_{1}\rho}\cos(\alpha-\beta)-\frac{\gamma_{1}Hu^{2}}{\rho
c}-\frac{\gamma_{1}Hv^{2}}{\rho c}+\frac{H^{2}\gamma_{1}}{\rho
c}\right)\frac{d\varphi_{\alpha-\beta}}{d\xi}\frac{d\xi}{ds}$ $\displaystyle($
$\displaystyle\textrm{using}\quad
u~{}\sin(\alpha-\beta)-v~{}\cos(\alpha-\beta)+c=0\quad\textrm{and}\quad
H=E_{c}+C^{2}/\gamma_{1})$ $\displaystyle\frac{d\Gamma^{1}_{S}}{ds}$
$\displaystyle=$
$\displaystyle\left(-\frac{c^{3}}{\gamma_{1}\rho}+\frac{\gamma_{1}}{\rho
c}(H-E_{c})^{2}\right)\frac{d\varphi_{\alpha-\beta}}{d\xi}\frac{d\xi}{ds}$
$\displaystyle($ $\displaystyle\textrm{using}\quad
u~{}\sin(\alpha-\beta)-v~{}\cos(\alpha-\beta)+c=0)$
$\displaystyle\frac{d\Gamma^{1}_{S}}{ds}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\frac{d\Gamma^{2}_{S}}{ds}$ $\displaystyle=$
$\displaystyle\frac{d\psi_{1}}{ds}+u\frac{d\psi_{2}}{ds}+v\frac{d\psi_{3}}{ds}+E_{c}\frac{d\psi_{4}}{ds}$
$\displaystyle\frac{d\Gamma^{2}_{S}}{ds}$ $\displaystyle=$
$\displaystyle\left(\frac{c}{\rho}(1+\frac{\gamma_{1}}{2}M^{2})+u(\frac{1}{\rho}(\sin(\alpha-\beta)-\gamma_{1}\displaystyle\frac{u}{c}))+v(\frac{1}{\rho}(-\cos(\alpha-\beta)-\gamma_{1}\displaystyle\frac{v}{c}))+E_{c}(\frac{\gamma_{1}}{\rho
c})\right)\frac{d\varphi_{\alpha-\beta}}{d\xi}\frac{d\xi}{ds}$
$\displaystyle\frac{d\Gamma^{2}_{S}}{ds}$ $\displaystyle=$
$\displaystyle\left(\frac{\gamma_{1}E_{c}}{\rho c}-\frac{\gamma_{1}u^{2}}{\rho
c}-\frac{\gamma_{1}v^{2}}{\rho c}+\frac{\gamma_{1}E_{c}}{\rho
c}\right)\frac{d\varphi_{\alpha-\beta}}{d\xi}\frac{d\xi}{ds}$ $\displaystyle($
$\displaystyle\textrm{using}\quad
u~{}\sin(\alpha-\beta)-v~{}\cos(\alpha-\beta)+c=0)$
$\displaystyle\frac{d\Gamma^{2}_{S}}{ds}$ $\displaystyle=$ $\displaystyle 0$

Besides

$\displaystyle\frac{d\Gamma_{C+}}{ds}$ $\displaystyle=$
$\displaystyle(u+vt^{+})\frac{d\psi_{1}}{ds}+(u^{2}+v^{2})\frac{d\psi_{2}}{ds}+(u^{2}+v^{2})t^{+}\frac{d\psi_{3}}{ds}+H(u+vt^{+})\frac{d\psi_{4}}{ds}$
$\displaystyle\frac{d\Gamma_{C+}}{ds}$ $\displaystyle=$
$\displaystyle((u+vt^{+})(\frac{c}{\rho}(1+\frac{\gamma_{1}}{2}M^{2}))+(u^{2}+v^{2})(\frac{1}{\rho}(\sin(\alpha-\beta)-\gamma_{1}\displaystyle\frac{u}{c}))$
$\displaystyle+(u^{2}+v^{2})t^{+}(\frac{1}{\rho}(-\cos(\alpha-\beta)-\gamma_{1}\displaystyle\frac{v}{c}))+H(u+vt^{+})(\frac{\gamma_{1}}{\rho
c}))\frac{d\varphi_{\alpha-\beta}}{d\xi}\frac{d\xi}{ds}$

It is easily checked that the first and fourth terms in the bracket are equal.
After lengthy calculations using the null eigenvalues properties and then the
formulas for the difference of $\cos$ and difference of $\sin$, it appears
that $d\Gamma_{C+}/ds=0$. For the sake of brevity, the detail of the
calculations is not shown here.

## Biblography

## References

  * [1] Jameson, A.  Aerodynamic design via control theory.  Journal of Scientific Computing 3(3), 233–260 (1988).
  * [2] Anderson, W. and Venkatakrishnan, V.  Adjoint design optimization on unstructured grid using a continuous adjoint formulation.  In AIAA Paper Series, Paper 1997-0643. (1997).
  * [3] Anderson, W., , and Venkatakrishnan, V.  Adjoint design optimization on unstructured grid using a continuous adjoint formulation.  Computers & Fluids 28, 443–480 (1999).
  * [4] Giles, M. and Pierce, N.  Adjoint equations in CFD: Duality, boundary conditions and solution behaviour.  In AIAA Paper Series, Paper 97-1850. (1997).
  * [5] Giles, M. and Pierce, N.  Analytic adjoint solutions for the quasi-one-dimensional Euler equations.  Journal of Fluid Mechanics 426, 327–345 (2001).
  * [6] Baeza, A., Castro, C., Palacios, F., and Zuazua, E.  2d Euler shape design on non-regular flows using adjoint Rankine-Hugoniot relations.  AIAA Journal 47(3), 552–562 (2009).
  * [7] Coquel, F. and Marmignon, C. and Rai, P. and Renac, F. Adjoint approximation of nonlinear hyperbolic systems with non-conservative products.  In XVII International Conference on Hyperbolic Problems: Theory, Numerics, Applications, Bressan, A., Lewicka, M., Wang, D., and Zheng, Y., editors, 385–392. AIMS, (2018).
  * [8] Lozano, C.  Singular and discontinuous solutions of the adjoint Euler equations.  AIAA Journal 56(11), 4437–4451 (2018).
  * [9] Peter, J., Renac, F., and Labbé, C.  Analysis of finite-volume discrete adjoint fields for two-dimensional compressible Euler flows.  Journal of Computational Physics 449, 110811 (2022).
  * [10] Lozano, C.  Watch your adjoints! lack of mesh convergence in inviscid adjoint solutions.  AIAA Journal 57(9), 3991–4006 (2019).
  * [11] Peter, J.  Contributions to discrete adjoint method in aerodynamics for shape optimization and goal-oriented mesh adaptation.  University of Nantes. Mémoire pour Habilitation à Diriger des Recherches, (2020).
  * [12] Ferri, A. Application of the method of characteristics to supersonic rotational flow.  Technical Report 841, NASA, (1946).
  * [13] Shapiro, A.H. The Dynamics and Thermodynamics of Compressible Fluid Flow.  The Ronald Press Company, (1954).
  * [14] Bonnet, A. and Luneau, J.  Aérodynamique. Théories de la dynamique des fluides. Cepadues, (1989).
  * [15] Délery, J.  Traité d’aérodynamique compressible. Volume 3.  Collection mécanique des fluides. Lavoisier – Hermés Science, (2008).
  * [16] Liepmann, H.W. and Roshko, A. Elements of Gasdynamics.  Dover Publications, (1956).
  * [17] Anderson, J.  Modern Compressible Flow (third edition).  McGraw-Hill series in Aernautical and Aerospace Engineering, (2003).
  * [18] Sears, W.R. (Under the editorial direction of).  General theory of high speed aerodynamics. Vol VI High speed aerodynamics and jet propulsion.  Princeton University Press, (1954).
  * [19] Venditti, D. and Darmofal, D.  Grid adaptation for functional outputs: Application to two-dimensional inviscid flows.  Journal of Computational Physics 176, 40–69 (2002).
  * [20] Todarello, G., Vonck, F., Bourasseau, S., Peter, J., and Désidéri, J.-A.  Finite-volume goal-oriented mesh-adaptation using functional derivative with respect to nodal coordinates.  Journal of Computational Physics 313, 799–819 (2016).
  * [21] Sartor, F. Unsteadiness in transonic shock-wave/boundary layer interactions: experimental investigation and global stability analysis.  PhD thesis, Université Aix-Marseille, (2014).
  * [22] Sartor, F. and Mettot, C. and Sipp, D. Stability, Receptivity, and Sensitivity Analyses of Buffeting Transonic Flow over a Profile.  AIAA Journal 53(7), 1980–1993 (2010).
  * [23] Hirsch, Ch. Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics (second edition).  Butterworth – Heineman. Elsevier, (2007).
  * [24] Jameson, A., Schmidt, W., and Turkel, E.  Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes.  In AIAA Paper Series, Paper 1981-1259. (1981).
  * [25] Peter, J. and Dwight, R.  Numerical sensitivity analysis for aerodynamic optimization: a survey of approaches.  Computers and Fluids 39, 373–391 (2010).
  * [26] Peter, J., Renac, F., Dumont, A., and Méheut, M.  Discrete adjoint method for shape optimization and mesh adaptation in the $els{A}$ code. status and challenges.  In Proceedings of 50th 3AF Symposium on Applied Aerodynamics, Toulouse, (2015).
  * [27] Mohammadi, B. and Pironneau, O.  Applied Shape Optimization for Fluids.  Oxford Univ. Press, May (2001).
  * [28] Alauzet, F. and Pironneau, O. Continuous and discrete adjoints to the Euler equations for fluids.  Int. J. Numer. Meth. Fluids 70, 135–157 (2012).
  * [29] Lozano, C. A note on the dual consistency of the discrete adjoint quasi-one dimensional Euler equations with cell-centred and cell-vertex central discretization.  Computers and Fluids 134-135, 51–60 (2016).
  * [30] Vassberg, J. and Jameson, A.  In pursuit of grid convergence for two-dimensional Euler solutions.  Journal of Aircraft 47(4), 1152–1166 (2010).