On Entropy-Stable Discretizations and the Entropy Adjoint

Abstract

For conservation laws that admit a convex entropy, Fidkowski and Roe (SIAM J Sci Comput 32(3):1261–1287, 2010) showed that the entropy variables satisfy an adjoint equation corresponding to an entropy-balance functional. In general, this relationship between the entropy variables and the entropy adjoint becomes an approximation when the conservation law is discretized. However, this work shows that the relationship is mimicked discretely for entropy-stable discretizations based on summation-by-parts (SBP) operators and two-point entropy-conservative flux functions; for these discretizations, the discrete entropy variables exactly satisfy the discrete adjoint equation for a particular discretization of the entropy-balance functional. A detailed proof of this result is presented for first-order conservation laws that are semi-discretized using SBP operators with diagonal norm and boundary operators. Subsequently, generalizations to second-order conservation laws, temporal discretizations, and non-diagonal SBP operators are described. The result is verified for steady inviscid and viscous flows.

Appendices

Alternative Form of the Discrete Residual

Here we show that $R_h({\mathbf {u}}_{h},{\mathbf {w}}_{h})$ defined by (24) is equivalent to the form in (25). First, substitute the strong form residual ${\mathbf {r}}_{i}(t)$ (see (16)) into the discrete residual:

$$\begin{aligned} R_h({\mathbf {u}}_{h},{\mathbf {w}}_{h})= & {} \int _{t=0}^T \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\\&\times \bigg [ \frac{d}{dt} {\mathbf {u}}_{i}+ \sum _{j=1}^{n} ({{\textsf {D}}}_{x})_{ij} 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) + \sum _{j=1}^{n} ({{\textsf {D}}}_{y})_{ij} 2 {\hat{{\mathbf {f}}}}_{y}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) \\&\quad + \sum _{j=1}^{n} {{\textsf {H}}}_{ii}^{-1} {{\textsf {M}}}_{ij}({\mathbf {u}}_{h}) {\mathbf {v}}_{j}- ({{\textsf {H}}}_{ii}^{-1} {{\textsf {E}}}_{ii}) \left[ {\mathbf {f}}_{n}({\mathbf {u}}_{i}) - {\mathbf {f}}_{\text {bc}}({\mathbf {u}}_{i}) \right] \bigg ] {{\textsf {H}}}_{ii} \, dt \\&+ \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}(0) ({\mathbf {u}}_{i}(0) - {\mathbf {u}}_{0,i}) {{\textsf {H}}}_{ii}. \end{aligned}$$

Next, we group the terms in a way that will be more convenient to handle individually:

$$\begin{aligned} R_h({\mathbf {u}}_{h},{\mathbf {w}}_{h})= & {} \int _{t=0}^T \left( \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\frac{d}{dt} {\mathbf {u}}_{i}\right) {{\textsf {H}}}_{ii} \, dt + \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}(0) ({\mathbf {u}}_{i}(0) - {\mathbf {u}}_{0,i}) {{\textsf {H}}}_{ii} \nonumber \\&\int _{t=0}^T \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\nonumber \\&\times \bigg [ \sum _{j=1}^{n} ({{\textsf {Q}}}_{x})_{ij} 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) + \sum _{j=1}^{n} ({{\textsf {Q}}}_{y})_{ij} 2 {\hat{{\mathbf {f}}}}_{y}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) - {{\textsf {E}}}_{ii} \left[ {\mathbf {f}}_{n}({\mathbf {u}}_{i}) - {\mathbf {f}}_{\text {bc}}({\mathbf {u}}_{i}) \right] \bigg ] \, dt \nonumber \\&+ \int _{t=0}^T \left( \sum _{i=1}^n \sum _{j=1}^{n} {\mathbf {w}}_{i}^{\mathsf {T}}{{\textsf {M}}}_{ij}({\mathbf {u}}_{h}) {\mathbf {v}}_{j}\right) \, dt, \end{aligned}$$

(54)

where we used ${{\textsf {H}}}_{ii} ({{\textsf {D}}}_{x})_{ij} = ({{\textsf {Q}}}_{x})_{ij}$ and ${{\textsf {H}}}_{ii} ({{\textsf {D}}}_{y})_{ij} = ({{\textsf {Q}}}_{y})_{ij}$. The last term on the right-hand side of (54) is already in the form presented in (25); we simply need to use the symmetry of ${{\textsf {M}}}_{ij}$ and swap the indices i and j. This leaves us to focus on the remaining terms.

Consider the first two terms on the right-hand-side of (54). Applying integration by parts to the time derivative term, we find

$$\begin{aligned}&\int _{t=0}^T \left( \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\frac{d}{dt} {\mathbf {u}}_{i}\right) {{\textsf {H}}}_{ii} \, dt + \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}(0) ({\mathbf {u}}_{i}(0) - {\mathbf {u}}_{0,i}) {{\textsf {H}}}_{ii} \nonumber \\&\quad = -\int _{t=0}^T \left( \sum _{i=1}^n {\mathbf {u}}_{i}^{\mathsf {T}}\frac{d}{dt}{\mathbf {w}}_{i}\right) {{\textsf {H}}}_{ii} \, dt + \sum _{i=1}^n \left[ \tilde{{\mathbf {u}}}_{i}^{\mathsf {T}}{\mathbf {w}}_{i}\right] _{t=0}^T {{\textsf {H}}}_{ii}, \end{aligned}$$

(55)

where, recall, $\tilde{{\mathbf {u}}}_{i}(0) \equiv {\mathbf {u}}_{0,i}$ and $\tilde{{\mathbf {u}}}_{i}(T) \equiv {\mathbf {u}}_{i}(T)$.

Finally, we address the spatial difference terms. For the moment, we ignore the temporal integral over the spatial difference terms in (54) and consider only the term containing the x direction entropy-conservative flux ${\hat{{\mathbf {f}}}}_{x}$:

$$\begin{aligned} \sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\sum _{j=1}^{n} ({{\textsf {Q}}}_{x})_{ij} 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j})&= \sum _{i=1}^n \sum _{j=1}^{n} ({{\textsf {Q}}}_{x})_{ij} {\mathbf {w}}_{i}^{\mathsf {T}}\left[ 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) - {\mathbf {f}}_x({\mathbf {u}}_{i})\right] \\&= \sum _{i=1}^n \sum _{j=1}^{n} ({{\textsf {Q}}}_{x})_{ij} {\mathbf {w}}_{i}^{\mathsf {T}}\tilde{{\mathbf {f}}}_{x,ij}, \end{aligned}$$

where we used the fact that the difference operator annihilates constants to introduce ${\mathbf {f}}_x({\mathbf {u}}_{i})$ on the first line; i.e. $\sum _{j=1}^n ({{\textsf {Q}}}_{x})_{ij} {\mathbf {f}}_x({\mathbf {u}}_{i}) = {\mathbf {0}}$. We have also introduced the numerical flux

$$\begin{aligned} \tilde{{\mathbf {f}}}_{x,ij} \equiv 2{\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) - {\mathbf {f}}_x({\mathbf {u}}_{i}). \end{aligned}$$

Continuing, we use the SBP property ${{\textsf {Q}}}_{x}= -{{\textsf {Q}}}_{x}^{\mathsf {T}}+ {{\textsf {E}}}_{x}$ to find

$$\begin{aligned}&\sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\sum _{j=1}^{n} ({{\textsf {Q}}}_{x})_{ij} 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) \nonumber \\&\quad = \sum _{i=1}^n \sum _{j=1}^{n} (-{{\textsf {Q}}}_{x}^{\mathsf {T}}+ {{\textsf {E}}}_{x})_{ij} {\mathbf {w}}_{i}^{\mathsf {T}}\tilde{{\mathbf {f}}}_{x,ij} \nonumber \\&\quad = -\sum _{i=1}^n \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ij}^{\mathsf {T}}({{\textsf {Q}}}_{x})_{ji} {\mathbf {w}}_{i}+ \sum _{i=1}^n {\mathbf {f}}_x^{\mathsf {T}}({\mathbf {u}}_{i}) {\mathbf {w}}_{i}\; n_{x,i} \; {{\textsf {E}}}_{ii} \qquad \text {(recall} \, {{\textsf {E}}}_{x}\, \text {is diagonal)} \nonumber \\&\quad = -\sum _{i=1}^n \left[ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {D}}}_{x})_{ij} {\mathbf {w}}_{j}\right] {{\textsf {H}}}_{ii} + \sum _{i=1}^n {\mathbf {f}}_x^{\mathsf {T}}({\mathbf {u}}_{i}) {\mathbf {w}}_{i}\; n_{x,i} \; {{\textsf {E}}}_{ii} \end{aligned}$$

(56)

where we used $\tilde{{\mathbf {f}}}_{x,ii} = {\mathbf {f}}_x({\mathbf {u}}_{i})$ in the boundary term, which follows from the consistency of the numerical flux: ${\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{i}) = {\mathbf {f}}_x({\mathbf {u}}_{i})$. Note, we also swapped the i and j indices in the double sum on the last line. This is significant, because the numerical flux $\tilde{{\mathbf {f}}}_{x,ij}$ is not symmetric in the indices.

Using the analogous identity to (56) for the y direction, we have the following:

$$\begin{aligned}&\sum _{i=1}^n {\mathbf {w}}_{i}^{\mathsf {T}}\bigg [ \sum _{j=1}^{n} ({{\textsf {Q}}}_{x})_{ij} 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) + \sum _{j=1}^{n} ({{\textsf {Q}}}_{y})_{ij} 2 {\hat{{\mathbf {f}}}}_{y}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) - {{\textsf {E}}}_{ii} \left[ {\mathbf {f}}_{n}({\mathbf {u}}_{i}) - {\mathbf {f}}_{\text {bc}}({\mathbf {u}}_{i}) \right] \bigg ] \nonumber \\&\quad = -\sum _{i=1}^n \left[ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {D}}}_{x})_{ij} {\mathbf {w}}_{j}+ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{y,ji}^{\mathsf {T}}({{\textsf {D}}}_{y})_{ij} {\mathbf {w}}_{j}\right] {{\textsf {H}}}_{ii} + {{\textsf {E}}}_{ii} {\mathbf {f}}_{\text {bc}}^{\mathsf {T}}({\mathbf {u}}_{i}) {\mathbf {w}}_{i}\end{aligned}$$

(57)

Substituting (55) and (57) into (54) produces the desired result, i.e. (25).

Derivative of the Spatial Difference Terms in the Discrete Residual

Lemma 3

Let $\tilde{{\mathbf {f}}}_{x,ji}$ and $\tilde{{\mathbf {f}}}_{y,ji}$ be the numerical fluxes defined by (26) and let ${{\textsf {D}}}_{x}= {{\textsf {H}}}^{-1} {{\textsf {Q}}}_{x}$ and ${{\textsf {D}}}_{y}= {{\textsf {H}}}^{-1} {{\textsf {Q}}}_{y}$ be diagonal-norm SBP operators in the x and y directions, respectively. Then,

$$\begin{aligned} \begin{aligned} \frac{\partial }{\partial {\mathbf {u}}_{k}} \sum _{i=1}^n \bigg [ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {D}}}_{x})_{ij} {\mathbf {w}}_{j}\bigg ] {{\textsf {H}}}_{ii}&= \bigg [ \sum _{j=1}^n 2 {\hat{{{\textsf {A}}}}}_{x,kj}^{\mathsf {T}}({{\textsf {D}}}_{x})_{kj} ({\mathbf {w}}_{j}- {\mathbf {w}}_{k}) \bigg ] {{\textsf {H}}}_{kk}, \\ \frac{\partial }{\partial {\mathbf {u}}_{k}} \sum _{i=1}^n \bigg [ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{y,ji}^{\mathsf {T}}({{\textsf {D}}}_{y})_{ij} {\mathbf {w}}_{j}\bigg ] {{\textsf {H}}}_{ii}&= \bigg [ \sum _{j=1}^n 2 {\hat{{{\textsf {A}}}}}_{y,kj}^{\mathsf {T}}({{\textsf {Q}}}_{y})_{kj} ({\mathbf {w}}_{j}- {\mathbf {w}}_{k}) \bigg ] {{\textsf {H}}}_{kk}. \end{aligned} \end{aligned}$$

(58)

Proof

The two identities in (58) are structurally identical, so it is sufficient to prove the first one. Differentiating with respect to ${\mathbf {u}}_{k}$, we find

$$\begin{aligned} \frac{\partial }{\partial {\mathbf {u}}_{k}} \sum _{i=1}^n \bigg [ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {D}}}_{x})_{ij} {\mathbf {w}}_{j}\bigg ] {{\textsf {H}}}_{ii}&= \sum _{i=1}^n \sum _{j=1}^n \left[ \frac{\partial }{\partial {\mathbf {u}}_{k}} 2 {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{i},{\mathbf {u}}_{j}) - \frac{\partial }{\partial {\mathbf {u}}_{k}} {\mathbf {f}}_x({\mathbf {u}}_{j})\right] ^T ({{\textsf {Q}}}_{x})_{ij} {\mathbf {w}}_{j}\\&= \sum _{i=1}^n \sum _{j=1}^n \left[ 2 {\hat{{{\textsf {A}}}}}_{x,ij} \delta _{ik} + 2{\hat{{{\textsf {A}}}}}_{x,ji} \delta _{jk} - {{\textsf {A}}}_{x}({\mathbf {u}}_{j}) \delta _{jk} \right] ^T ({{\textsf {Q}}}_{x})_{ij} {\mathbf {w}}_{j}, \end{aligned}$$

where $\delta _{ik}$ and $\delta _{jk}$ denote Kronecker deltas. Continuing,

$$\begin{aligned} \frac{\partial }{\partial {\mathbf {u}}_{k}} \sum _{i=1}^n \bigg [ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {D}}}_{x})_{ij} {\mathbf {w}}_{j}\bigg ] {{\textsf {H}}}_{ii}&= \sum _{j=1}^n 2 {\hat{{{\textsf {A}}}}}_{x,kj}^{\mathsf {T}}({{\textsf {Q}}}_{x})_{kj} {\mathbf {w}}_{j}+ \sum _{i=1}^n \left[ 2 {\hat{{{\textsf {A}}}}}_{x,ki}^{\mathsf {T}}- {{\textsf {A}}}_x^{\mathsf {T}}({\mathbf {u}}_{k}) \right] ({{\textsf {Q}}}_{x})_{ik} {\mathbf {w}}_k \\&=\sum _{j=1}^n 2 {\hat{{{\textsf {A}}}}}_{x,kj}^{\mathsf {T}}({{\textsf {Q}}}_{x})_{kj} {\mathbf {w}}_{j}\\&\quad + \sum _{j=1}^n \left[ 2 {\hat{{{\textsf {A}}}}}_{x,kj}^{\mathsf {T}}- {{\textsf {A}}}_x^{\mathsf {T}}({\mathbf {u}}_{k}) \right] (-{{\textsf {Q}}}_{x}^{\mathsf {T}}+ {{\textsf {E}}}_{x})_{jk} {\mathbf {w}}_k, \end{aligned}$$

where we have introduced the SBP property ${{\textsf {Q}}}_{x}= -{{\textsf {Q}}}_{x}^{\mathsf {T}}+ {{\textsf {E}}}_{x}$, and we replaced the dummy index i with j in the second sum. Rearranging terms we find

$$\begin{aligned} \frac{\partial }{\partial {\mathbf {u}}_{k}} \sum _{i=1}^n \bigg [ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {D}}}_{x})_{ij} {\mathbf {w}}_{j}\bigg ] {{\textsf {H}}}_{ii}&= \sum _{j=1}^n 2 {\hat{{{\textsf {A}}}}}_{x,kj}^{\mathsf {T}}({{\textsf {Q}}}_{x})_{kj} ({\mathbf {w}}_{j}- {\mathbf {w}}_{k}) \\&\quad + {{\textsf {A}}}_x^{\mathsf {T}}({\mathbf {u}}_{k}) {\mathbf {w}}_k \sum _{j=1}^n ({{\textsf {Q}}}_{x})_{kj} + \left[ 2 {\hat{{{\textsf {A}}}}}_{x,kk}^{\mathsf {T}}- {{\textsf {A}}}_x^{\mathsf {T}}({\mathbf {u}}_{k})\right] {\mathbf {w}}_k \; n_{x,k} {{\textsf {E}}}_{kk}. \end{aligned}$$

The two last terms on the right-hand side above are zero, because $\sum _{j=1}^n ({{\textsf {Q}}}_{x})_{kj} = 0$ and because, as shown in, for example, the proof of Theorem 1 of [6],

$$\begin{aligned} {{\textsf {A}}}_x({\mathbf {u}}_{k}) = \frac{\partial }{\partial {\mathbf {u}}_{k}} {\mathbf {f}}_x({\mathbf {u}}_{k}) = \frac{\partial }{\partial {\mathbf {u}}_{k}} {\hat{{\mathbf {f}}}}_{x}({\mathbf {u}}_{k},{\mathbf {u}}_{k}) = 2 {\hat{{{\textsf {A}}}}}_{x,kk}. \end{aligned}$$

Therefore, we are left with

$$\begin{aligned} \frac{\partial }{\partial {\mathbf {u}}_{k}} \sum _{i=1}^n \bigg [ \sum _{j=1}^n \tilde{{\mathbf {f}}}_{x,ji}^{\mathsf {T}}({{\textsf {Q}}}_{x})_{ij} {\mathbf {w}}_{j}\bigg ] {{\textsf {H}}}_{ii} = \sum _{j=1}^n 2 {\hat{{{\textsf {A}}}}}_{x,kj}^{\mathsf {T}}({{\textsf {Q}}}_{x})_{kj} ({\mathbf {w}}_{j}- {\mathbf {w}}_{k}). \end{aligned}$$

(59)

The result follows after introducing ${{\textsf {H}}}_{kk}{{\textsf {H}}}_{kk}^{-1}$ on the right-hand side. $\square$

Discrete Entropy-Adjoint for Second-Order Conservation Laws

As we did for the discrete adjoint analysis of first-order conservation laws, we begin by defining the weak formulation. For the discrete second-order conservation law we have

$$\begin{aligned}&R_h({\mathbf {u}}_{h},{\mathbf {w}}_{h}) \equiv \int _{t=0}^{T} \sum _{i=1}^{n} {\mathbf {w}}_i^{\mathsf {T}}{\mathbf {r}}_i(t) {{\textsf {H}}}_{ii} \, dt + \sum _{i=1}^{n} {\mathbf {w}}_{i}^{\mathsf {T}}(0) ({\mathbf {u}}_{i}(0) - {\mathbf {u}}_{0,i}) {{\textsf {H}}}_{ii} = {\mathbf {0}}, \nonumber \\&\qquad \forall \; {\mathbf {w}}_{h}\in {\mathbb {R}}^{ns}, \end{aligned}$$

(60)

where ${\mathbf {r}}_{i}(t)$ is defined by (40). The temporal and inviscid terms that make up $R_h({\mathbf {u}}_{h},{\mathbf {w}}_{h})$ have previously been shown to be equivalent to (25). Therefore, here we focus on the viscous terms.

Consider the ${{\textsf {H}}}$ inner product between the test function ${\mathbf {w}}_{h}$ and the discrete viscous terms defined by (41). This product is equal to (recall the SBP operator definitions ${{\textsf {D}}}_{x}= {{\textsf {H}}}^{-1} {{\textsf {Q}}}_{x}$ and ${{\textsf {D}}}_{y}= {{\textsf {H}}}^{-1} {{\textsf {Q}}}_{y}$)

$$\begin{aligned} \sum _{i=1}^n {\mathbf {w}}_i^{\mathsf {T}}{\mathbf {r}}_i^V(t) {{\textsf {H}}}_{ii}= & {} -\sum _{i=1}^n \sum _{j=1}^n {\mathbf {w}}_i^{\mathsf {T}}\left[ ({{\textsf {Q}}}_{x})_{ij} {\mathbf {f}}_{x,j}^V({\mathbf {u}}_{h}) + ({{\textsf {Q}}}_{y})_{ij} {\mathbf {f}}_{y,j}^V({\mathbf {u}}_{h}) \right] \\&+ \sum _{i=1}^n {\mathbf {w}}_i^{\mathsf {T}}{{\textsf {E}}}_{ii} \left[ {\mathbf {f}}_{n,i}^V({\mathbf {u}}_{h}) - {\mathbf {f}}_{\text {bc},i}^V({\mathbf {u}}_{h}) \right] . \end{aligned}$$

In order to simplify the weak-form of the viscous terms, we use the SBP property and the definition of the viscous flux in the normal direction; see Eq. (44).

$$\begin{aligned} \sum _{j=1}^n {\mathbf {w}}_i^{\mathsf {T}}{\mathbf {r}}_i^V(t) {{\textsf {H}}}_{ii} = \sum _{i=1}^n \sum _{j=1}^n {\mathbf {w}}_i^{\mathsf {T}}\left[ ({{\textsf {Q}}}_{x})_{ji} {\mathbf {f}}_{x,j}^V({\mathbf {u}}_{h}) + ({{\textsf {Q}}}_{y})_{ji} {\mathbf {f}}_{y,j}^V({\mathbf {u}}_{h}) \right] - \sum _{i=1}^n {\mathbf {w}}_i^{\mathsf {T}}{{\textsf {E}}}_{ii} {\mathbf {f}}_{\text {bc},i}^V({\mathbf {u}}_{h}). \end{aligned}$$

Finally, we substitute the definitions of ${\mathbf {f}}_{x,j}^V({\mathbf {u}}_{h})$ and ${\mathbf {f}}_{y,j}^V({\mathbf {u}}_{h})$ from (42) and (43), respectively, to obtain

$$\begin{aligned} \sum _{j=1}^n {\mathbf {w}}_i^{\mathsf {T}}{\mathbf {r}}_i^V(t) {{\textsf {H}}}_{ii} = \sum _{i=1}^n \sum _{j=1}^n {\mathbf {w}}_i^{\mathsf {T}}{{\textsf {M}}}_{ij}^V {\mathbf {v}}_{j}- \sum _{i=1}^n {\mathbf {w}}_i^{\mathsf {T}}{\mathbf {f}}_{\text {bc},i}^V({\mathbf {u}}_{h}) {{\textsf {E}}}_{ii}, \end{aligned}$$

where the entries in the symmetric, positive semi-definite matrix ${{\textsf {M}}}^V$ are defined by

$$\begin{aligned} {{\textsf {M}}}_{ij}^V({\mathbf {u}}_{h})\equiv & {} \sum _{k=1}^n ({{\textsf {D}}}_{x})_{ki} {{\textsf {H}}}_{kk} \left[ ({{\textsf {C}}}_{xx})_k ({{\textsf {D}}}_{x})_{kj} + ({{\textsf {C}}}_{xy})_k ({{\textsf {D}}}_{y})_{kj} \right] \nonumber \\&+ \sum _{k=1}^n ({{\textsf {D}}}_{y})_{ki} {{\textsf {H}}}_{kk} \left[ ({{\textsf {C}}}_{yx})_k ({{\textsf {D}}}_{x})_{kj} + ({{\textsf {C}}}_{yy})_k ({{\textsf {D}}}_{y})_{kj} \right] . \end{aligned}$$

(61)

Thus, the semi-linear form in the weak form (60) can be written as

$$\begin{aligned} R_h({\mathbf {u}}_{h},{\mathbf {w}}_{h})= & {} (\text {terms on RHS of} \, (25)) - \int _{t=0}^T \sum _{i=1}^n {\mathbf {w}}_i^{\mathsf {T}}{\mathbf {f}}_{\text {bc},i}^V({\mathbf {u}}_{h}) {{\textsf {E}}}_{ii} \, dt \nonumber \\&+ \int _{t=0}^T \sum _{i=1}^n \sum _{j=1}^n {\mathbf {w}}_i^{\mathsf {T}}{{\textsf {M}}}_{ij}^V {\mathbf {v}}_{j}\, dt. \end{aligned}$$

(62)

The discrete adjoint analysis for the second-order conservation law is analogous to the analysis presented earlier in Sect. 3.4. Differentiating $R_h({\mathbf {u}}_{h}, {\mathbf {w}}_{h})$ with respect to ${\mathbf {u}}_{h}$ and replacing the test function with the entropy variables, we find (recall that ${{\textsf {M}}}_{ij}^V$ is a function of ${\mathbf {u}}_{h}$, in general)

$$\begin{aligned} R_h'[{\mathbf {u}}_{h}](\delta {\mathbf {u}}_{h},{\mathbf {v}}_{h})= & {} (\text {terms on RHS of} \, (27)) - \int _{t=0}^T \sum _{i=1}^n {\mathbf {v}}_{i}^{\mathsf {T}}\left( \frac{\partial {\mathbf {f}}_{\text {bc},i}^V}{\partial {\mathbf {u}}_{h}} \right) \delta {\mathbf {u}}_{h}{{\textsf {E}}}_{ii} \, d t \nonumber \\&+ \int _{t=0}^T \sum _{i=1}^n \sum _{j=1}^n \left[ {\mathbf {v}}_{i}^{\mathsf {T}}{{\textsf {M}}}_{ij}^V {{\textsf {A}}}_{0,j}^{-1} \delta {\mathbf {u}}_{j}+ \left( {\mathbf {v}}_{i}^{\mathsf {T}}\frac{\partial {{\textsf {M}}}_{ij}^V}{\partial {\mathbf {u}}_{h}} {\mathbf {v}}_{j}\right) \delta {\mathbf {u}}_{h}\right] \, dt, \end{aligned}$$

(63)

where ${{\textsf {A}}}_{0,j}^{-1} = \partial {\mathbf {v}}_{j}/\partial {\mathbf {u}}_{j}$ is symmetric. Next, we add and subtract terms from $R_h'[{\mathbf {u}}_{h}](\delta {\mathbf {u}}_{h}, {\mathbf {v}}_{h})$ in order to obtain the following expression.

$$\begin{aligned} R_h'[{\mathbf {u}}_{h}](\delta {\mathbf {u}}_{h},{\mathbf {v}}_{h})= & {} - \int _{t=0}^T \sum _{i=1}^n \delta {\mathbf {u}}_{i}^{\mathsf {T}}{{\textsf {A}}}_{0,i}^{-1} {\mathbf {r}}_{i}(t) {{\textsf {H}}}_{ii} \, dt + (\text {terms on RHS of} \, (29)) \nonumber \\&+ \underbrace{\int _{t=0}^T \sum _{i=1}^n \delta {\mathbf {u}}_{i}^{\mathsf {T}}{{\textsf {A}}}_{0,i}^{-1} \left[ \sum _{j=1}^n {{\textsf {H}}}_{ii}^{-1} {{\textsf {M}}}_{ij}^V {\mathbf {v}}_{j}- ({{\textsf {H}}}_{ii}^{-1} {{\textsf {E}}}_{ii}) {\mathbf {f}}_{\text {bc},i}^V({\mathbf {u}}_{h}) \right] {{\textsf {H}}}_{ii} \, dt}_{\displaystyle \text {added to balance terms introduced for} \, {\mathbf {r}}_{i}(t)} \nonumber \\&- \int _{t=0}^T \sum _{i=1}^n {\mathbf {v}}_{i}^{\mathsf {T}}\left( \frac{\partial {\mathbf {f}}_{\text {bc},i}^V}{\partial {\mathbf {u}}_{h}} \right) \delta {\mathbf {u}}_{h}\, d t \nonumber \\&+ \int _{t=0}^T \sum _{i=1}^n \sum _{j=1}^n \left[ {\mathbf {v}}_{i}^{\mathsf {T}}{{\textsf {M}}}_{ij}^V {{\textsf {A}}}_{0,j}^{-1} \delta {\mathbf {u}}_{j}+ \left( {\mathbf {v}}_{i}^{\mathsf {T}}\frac{\partial {{\textsf {M}}}_{ij}^V}{\partial {\mathbf {u}}_{h}} {\mathbf {v}}_{j}\right) \delta {\mathbf {u}}_{h}\right] \, dt. \end{aligned}$$

(64)

Grouping boundary and volume terms—and recognizing that ${\mathbf {r}}_i(t) = {\mathbf {0}}$ in light of (40)—we arrive at

$$\begin{aligned}&R_h'[{\mathbf {u}}_{h}](\delta {\mathbf {u}}_{h},{\mathbf {v}}_{h}) \nonumber \\&\quad = (\text {terms on RHS of} \, (29)) \nonumber \\&\qquad - \int _{t=0}^T \sum _{i=1}^n \left[ \delta {\mathbf {u}}_{i}^{\mathsf {T}}{{\textsf {A}}}_{0,i}^{-1} {\mathbf {f}}_{\text {bc},i}^V({\mathbf {u}}_{h}) + {\mathbf {v}}_{i}^{\mathsf {T}}\left( \frac{\partial {\mathbf {f}}_{\text {bc},i}^V}{\partial {\mathbf {u}}_{h}} \right) \delta {\mathbf {u}}_{h}\right] {{\textsf {E}}}_{ii} \, dt \nonumber \\&\qquad + \int _{t=0}^T \sum _{i=1}^n \sum _{j=1}^n \left[ {\mathbf {v}}_{i}^{\mathsf {T}}{{\textsf {M}}}_{ij}^V {{\textsf {A}}}_{0,j}^{-1} \delta {\mathbf {u}}_{j}+ \delta {\mathbf {u}}_{i}^T {{\textsf {A}}}_{0,i}^{-1} {{\textsf {M}}}_{ij}^V {\mathbf {v}}_{j}+ \left( {\mathbf {v}}_{i}^{\mathsf {T}}\frac{\partial {{\textsf {M}}}_{ij}^V}{\partial {\mathbf {u}}_{h}} {\mathbf {v}}_{j}\right) \delta {\mathbf {u}}_{h}\right] \, dt. \end{aligned}$$

(65)

Consequently, if we define the discretized entropy functional as

then $R_h'[{\mathbf {u}}_{h}](\delta {\mathbf {u}}_{h}, {\mathbf {v}}_{h}) = J_h'[{\mathbf {u}}_{h}](\delta {\mathbf {u}}_{h})$, as desired.