Discrete adjoints for one-dimensional laminar Navier

阅读范围：教学主线与研究性选读

Reading scope: core teaching and optional research

本项目以离散伴随教学为主。第 1–15、17–18 节构成从残差、tangent、adjoint 到梯度验证与线性求解的主线；第 16 节先了解“离散导数正确”与“连续极限正确”的区别即可。附录 A 的独立连续参考和附录 C 的直接法是进阶教学；附录 B 逐面核对第 8、9 节的局部求导。

The teaching core is Sections 1–15 and 17–18: residuals, tangents, adjoints, gradient checks and linear solves. In Section 16, first distinguish correct discrete derivatives from a correct continuum limit. The independent continuous reference in Appendix A and the direct solve in Appendix C are advanced teaching; Appendix B checks the local derivatives of Sections 8 and 9 face by face.

以下标记只表示超出首次学习的研究性分析范围，不表示首次提出、学术空白或已经证明的新结论。可以跳过这些分析而完整学习离散伴随。

The labels identify analysis beyond a first course, not novelty, a literature gap or a newly proved result. These parts can be skipped without missing the discrete-adjoint method.

离散梯度的连续极限Continuum limits of discrete gradients

第 1 节Section 1

方程：一维定常层流 Navier–Stokes 方程组

The equations: steady one-dimensional laminar Navier–Stokes

先明确这个算例求什么

What this model actually solves

本页求解等截面、平面一维、可压缩层流的稳态问题。未知量只随轴向位置变化，横向速度为零；黏性来自分子应力，传热来自 Fourier 热传导。它不解析管壁边界层，不含湍流模型，也不是变截面喷管的黏性修正。教学重点是：同一段残差程序如何产生原始解、精确 tangent 乘积和精确转置乘积。

This page solves steady, compressible, plane one-dimensional laminar flow at constant area. Unknowns vary only axially and transverse velocity is zero. Molecular stress and Fourier conduction supply viscous transport. It resolves neither duct-wall boundary layers nor turbulence and is not a viscous correction to a variable-area nozzle. The teaching focus is how one residual program supplies the primal solution, exact tangent products and exact transpose products.

从守恒律到定常方程组

From conservation laws to the steady equations

一维可压缩流有三条守恒律：质量、动量和总能量。用守恒变量 $U=(\rho,\rho u,\rho E)^T$ 表示时，任一区间内这三个量的变化只来自区间两端的通量和区间内的源项。写成微分形式，就是一维 Navier–Stokes 方程组：

One-dimensional compressible flow has three conservation laws: mass, momentum and total energy. In the conservative variables $U=(\rho,\rho u,\rho E)^T$ , the amount of each in an interval changes only through the fluxes at its two ends and the source inside it. In differential form these are the one-dimensional Navier–Stokes equations:

U_t+\bigl(F^c(U)-F^v\bigr)_x=S(x,D),\qquad F^c=\begin{bmatrix}\rho u\\\rho u^2+p\\u(\rho E+p)\end{bmatrix},\qquad F^v=\begin{bmatrix}0\\\tau\\u\tau-q\end{bmatrix}.

这里 $U(x,t)$ 是连续的守恒场， $E$ 是单位质量的总能量；第 2 节起的 $U_i$ 、 $\mathbf U_h$ 是它的离散版本，附录 A 需要区分时把连续场写作 $U^{\rm c}$ 。对流通量 $F^c$ 与 Euler 方程相同：流体携带质量、动量和能量，压力做功。新增的是黏性通量 $F^v$ ：黏性正应力 $\tau$ 传递动量，黏性功 $u\tau$ 与热流 $q$ 传递能量；总通量是对流减黏性。状态和物性由下面的关系闭合：

Here $U(x,t)$ is the continuous conservative field and $E$ the total energy per unit mass; the $U_i$ and $\mathbf U_h$ of Section 2 onwards are its discrete versions, and Appendix A writes the continuous field $U^{\rm c}$ where the two must be told apart. The convective flux $F^c$ is that of the Euler equations: the fluid carries mass, momentum and energy, and the pressure does work. New is the viscous flux $F^v$ : the viscous normal stress $\tau$ transfers momentum, and the viscous work $u\tau$ and the heat flux $q$ transfer energy; the total flux is convection minus viscosity. The state and the transport properties are closed by

\begin{gathered}E=c_vT+\tfrac12u^2,\quad p=\rho T,\quad c=\sqrt{\gamma T},\qquad\tau=\frac43\mu u_x,\quad q=-\kappa T_x,\\\mu=\mu_{\rm ref}\ \ \text{(constant)}\quad\text{or}\quad\mu=\mu_{\rm ref}\frac{1.4\,T^{3/2}}{T+0.4}\ \ \text{(Sutherland)},\qquad\kappa=\frac{c_p}{Pr}\mu.\end{gathered}

采用无量纲气体常数 1， $\gamma=1.4$ 、 $c_v=2.5$ 、 $c_p=3.5$ 。应力系数 $4/3$ 来自三维 Newton 应力在平面一维运动下的 Stokes 假设。质量分量没有分子扩散。能量黏性通量同时包含黏性功和导热，漏掉任意一项都会改变原始问题和它的导数。Sutherland 式中的 $1.4=1+S_\mu$ ， $S_\mu=0.4$ 是无量纲 Sutherland 常数，使 $\mu(1)=\mu_{\rm ref}$ ；它不是比热比 $\gamma$ 。输入的 $\mu_{\rm ref}$ 已按参考密度、长度和速度无量纲化，不再额外乘 $1/Re$ 。

The nondimensional gas constant is one, with $\gamma=1.4$ , $c_v=2.5$ , $c_p=3.5$ . The factor $4/3$ follows from three-dimensional Newtonian stress under Stokes' hypothesis restricted to plane motion. There is no molecular mass diffusion. The viscous energy flux contains both viscous work and conduction; omitting either changes the primal problem and its derivative. In the Sutherland law, $1.4=1+S_\mu$ with the nondimensional Sutherland constant $S_\mu=0.4$ , so that $\mu(1)=\mu_{\rm ref}$ ; it is not the ratio of specific heats $\gamma$ . The input $\mu_{\rm ref}$ is already scaled by reference density, length and velocity; no additional $1/Re$ is applied.

本页只研究稳态解。令 $U_t=0$ ，就得到要离散的方程组和它的边界条件（页顶默认的加热算例）：

This page studies steady solutions only. Setting $U_t=0$ gives the equations to be discretised, with their boundary conditions (the heated case, the default at the top of the page):

\begin{gathered}\frac{d}{dx}\bigl(F^c(U)-F^v\bigr)=S(x,D),\qquad x\in[0,1],\\\rho=1,\quad u=0.35,\quad T=T_{\rm in}\quad\text{at }x=0;\qquad p-\tau=p_b,\quad q=0\quad\text{at }x=1.\end{gathered}

解是光滑的，处处 $u>0$ 且保持亚声速，没有激波，所以第 3 节的二阶格式使用不加 limiter 的线性重构。

The solution is smooth, with $u>0$ everywhere, and stays subsonic, without shocks, so the second-order scheme of Section 3 uses linear reconstruction without a limiter.

图 1把这个问题画在一张图上。上面是流道：热源 $S_E(x)$ 集中在区间中部，入口给定 3 个量，出口给定 2 个量。下面四幅小图是默认参数下的解：曲线是附录 A 的独立连续参考解（256 步 RK4 的射击法），圆点是本页算出的离散解（ $N=12$ 、二阶，节点值；第 2 节起讲怎样算）。四个量相对入口的变化都不到 1%，所以每幅图的纵轴只画很窄的一段。圆点与曲线之差是 $N=12$ 时的离散误差，四个量中最大约 $2.9\times10^{-4}$ ，网格加密一倍后降到约 $1.2\times10^{-4}$ 。最下面一行是六个设计分量的当前值。

Figure 1 draws the whole problem. On top is the duct: the heat source $S_E(x)$ is concentrated in the middle of the interval, 3 quantities are prescribed at the inlet and 2 at the outlet. The four small panels below show the solution at the default parameters: the curves are the independent continuous reference of Appendix A (shooting with 256 RK4 steps), the dots the discrete solution this page computes ( $N=12$ , second order, nodal values; Section 2 onwards shows how). Each quantity changes by less than 1% from its inlet value, so each vertical axis spans a narrow range. The dots differ from the curves by the discretisation error at $N=12$ , at most about $2.9\times10^{-4}$ over the four quantities, and about $1.2\times10^{-4}$ on a mesh twice as fine. The bottom line gives the current values of the six design components.

图 1 · 问题全貌（页顶默认参数：加热算例，N = 12，二阶，节点值）。上：流道与热源 S_E(x) = a₀ sin²(πx) + a₁ sin²(2πx)，入口给定 ρ、u、T 共 3 个条件，出口给定 p − τ 与 q 共 2 个条件。下：温度、密度、速度与压力；曲线是附录 A 的独立连续参考解（256 步 RK4），圆点是本页的离散解，与曲线最多相差约 2.9×10⁻⁴。四个量相对入口的变化都小于 1%，所以各自的纵轴只画很窄的一段。底部：六个设计分量的默认值。

Figure 1 · The problem at a glance (the defaults at the top of the page: heated case, N = 12, second order, nodal values). Top: the duct and the heat source S_E(x) = a₀ sin²(πx) + a₁ sin²(2πx); the inlet prescribes 3 conditions (ρ, u, T), the outlet 2 (p − τ and q). Below: temperature, density, velocity and pressure; the curves are the independent continuous reference of Appendix A (256 RK4 steps), the dots this page's discrete solution, at most about 2.9×10⁻⁴ away from them. Each changes by less than 1% from its inlet value, so each vertical axis spans a narrow range. Bottom: the default values of the six design components.

波速、特征与两端的条件个数

Wave speeds, characteristics and the conditions at each end

先看没有黏性和导热的部分，即 Euler 方程。对流通量的 Jacobian $\partial F^c/\partial U$ 有三个特征值 $u-c$ 、 $u$ 、 $u+c$ ，其中 $c=\sqrt{\gamma T}$ 是声速：小扰动沿这三个速度传播，它们叫作特征速度。本页的流动处处满足 $0<u<c$ （入口 $\mathrm{Ma}=u/c\approx0.30$ ），所以 $u-c<0<u<u+c$ ：入口处三者依次为 $-0.833$ 、 $0.350$ 、 $1.533$ 。在入口， $u$ 与 $u+c$ 两族特征进入区间， $u-c$ 一族离开；在出口只有 $u-c$ 一族进入（图 2(a)）。每一族进入的特征需要一个边界条件，所以 Euler 方程在入口给 2 个条件、在出口给 1 个，拟一维 Euler 页正是这样（入口给总压和总温，出口给静压）。

Start with the part without viscosity and conduction, the Euler equations. The Jacobian $\partial F^c/\partial U$ of the convective flux has three eigenvalues, $u-c$ , $u$ and $u+c$ , where $c=\sqrt{\gamma T}$ is the speed of sound: small disturbances travel with these three speeds, the characteristic speeds. The flow on this page has $0<u<c$ everywhere (at the inlet $\mathrm{Ma}=u/c\approx0.30$ ), so $u-c<0<u<u+c$ : at the inlet they are $-0.833$ , $0.350$ and $1.533$ . At the inlet two families, $u$ and $u+c$ , enter the interval and $u-c$ leaves; at the outlet only $u-c$ enters (Figure 2(a)). Each entering family needs one boundary condition, so the Euler equations take 2 conditions at the inlet and 1 at the outlet, as on the quasi-1-D Euler page (total pressure and total temperature at the inlet, static pressure at the outlet).

黏性和导热改变了这个计数。 $\tau=\frac43\mu u_x$ 与 $q=-\kappa T_x$ 把 $u$ 和 $T$ 的二阶导数带进动量方程和能量方程，这两个方程于是有了抛物型（扩散）的部分：扰动不只沿特征传播，还同时向上游和下游扩散。定常时，动量方程对 $u$ 、能量方程对 $T$ 都是二阶的，各在两端要一个条件；质量方程没有扩散项，仍是一阶，由 $u>0$ 从入口输运，只在入口要一个条件。合起来，定常方程组的阶数是 $1+2+2=5$ ：入口 3 个条件，出口 2 个（图 2(b)）。本页入口给定 $\rho$ 、 $u$ 、 $T$ ，出口给定法向总应力 $p-\tau$ 与热流 $q$ ：动量方程的两个条件是入口的 $u$ 与出口的 $p-\tau$ ，能量方程的是入口的 $T$ 与出口的 $q$ ，质量方程的是入口的 $\rho$ 。

Viscosity and conduction change this count. $\tau=\frac43\mu u_x$ and $q=-\kappa T_x$ put second derivatives of $u$ and $T$ into the momentum and energy equations, which therefore gain a parabolic (diffusive) part: disturbances no longer travel only along characteristics but also diffuse upstream and downstream at the same time. In the steady equations the momentum equation is second order in $u$ and the energy equation second order in $T$ , so each takes one condition at each end; the mass equation has no diffusion, stays first order and is carried from the inlet by $u>0$ , so it takes one condition, at the inlet. In total the steady system has order $1+2+2=5$ : 3 conditions at the inlet and 2 at the outlet (Figure 2(b)). This page prescribes $\rho$ , $u$ and $T$ at the inlet and the total normal stress $p-\tau$ and the heat flux $q$ at the outlet: the momentum equation's two conditions are $u$ at the inlet and $p-\tau$ at the outlet, the energy equation's are $T$ at the inlet and $q$ at the outlet, and the mass equation's is $\rho$ at the inlet.

图 2 · 两端各要几个条件。(a) Euler 部分：本页离散解在入口面和出口面的特征速度 u − c、u、u + c，画成 x–t 平面上从边界点出发的射线。实线进入区间，需要一个边界条件；虚线离开区间，由内部决定。入口有 2 族进入、出口有 1 族进入。(b) 加上黏性与导热后，动量方程对 u、能量方程对 T 成为二阶，各在两端要一个条件；质量方程仍是一阶，只在入口要一个。合计入口 3 个、出口 2 个，格点页正好替换 3 + 2 行残差。

Figure 2 · How many conditions each end takes. (a) The Euler part: the characteristic speeds u − c, u and u + c of this page's discrete solution at the inlet and outlet faces, drawn as rays from a boundary point in the x–t plane. Solid rays enter the interval and need a boundary condition; dashed rays leave it and are set by the interior. 2 families enter at the inlet and 1 at the outlet. (b) With viscosity and conduction the momentum equation becomes second order in u and the energy equation second order in T, so each takes one condition at each end; the mass equation stays first order and takes one, at the inlet. In total 3 at the inlet and 2 at the outlet: the node page replaces exactly 3 + 2 residual rows.

特征速度 $u\pm c$ 并没有失去作用：它们仍然描述对流部分，第 3 节数值耗散的速度 $\alpha_f$ 就是面两侧 $u+c$ 的平均。默认参数下 $Re=\rho uL/\mu_{\rm ref}\approx2.3$ 、 $Pe=Re\cdot Pr\approx1.7$ ，都是 1 的量级：黏性和导热在整个区间都起作用，而不只在两端的薄层里。数出条件的个数只说明每端要给几个量，它本身不证明问题适定；第 4 节把这些条件写成实际的边界映射。

The characteristic speeds $u\pm c$ keep their role: they still describe the convective part, and the dissipation speed $\alpha_f$ of Section 3 is the average of $u+c$ on the two sides of a face. With the default parameters $Re=\rho uL/\mu_{\rm ref}\approx2.3$ and $Pe=Re\cdot Pr\approx1.7$ are of order one: viscosity and conduction act over the whole interval, not only in thin layers at the ends. Counting conditions only says how many quantities each end needs; it does not by itself prove the problem well-posed, and Section 4 turns these conditions into the implemented boundary maps.

本页的问题：数据、设计变量与基准解

The problem on this page: data, design variables and the baseline solution

页顶的默认算例是一段受热的等截面流道，区间为 $[0,1]$ ：入口 $x=0$ 给定 $\rho=1$ 、 $u=0.35$ 、 $T=T_{\rm in}$ ，出口 $x=1$ 给定 $p-\tau=p_b$ 与 $q=0$ 。热源只进入能量方程：

The default case at the top of the page is a heated duct of constant area on $[0,1]$ : the inlet $x=0$ prescribes $\rho=1$ , $u=0.35$ and $T=T_{\rm in}$ , the outlet $x=1$ prescribes $p-\tau=p_b$ and $q=0$ . The heat source enters the energy equation only:

S(x,D)=\begin{bmatrix}0\\0\\S_E(x,D)\end{bmatrix},\qquad S_E(x,D)=a_0\sin^2(\pi x)+a_1\sin^2(2\pi x),\qquad\int_0^1S_E\,dx=\frac{a_0+a_1}{2}.

设计变量是六维向量 $D=(a_0,a_1,\mu_{\rm ref},Pr,T_{\rm in},b_{\rm out})$ ：两个热源幅度，两个物性参数（参考黏度与 Prandtl 数），两个边界数据——入口温度 $T_{\rm in}=D_4$ 与出口数据 $b_{\rm out}=D_5$ （本算例是 $p_b$ ）。按伴随优化文献的通常写法，设计变量记作 $D$ 。第 5 节把 $D$ 写成离散残差的第二个自变量。默认值是 $D=(0.03,\,0,\,0.15,\,0.72,\,1,\,1)$ （图 1底部），加入的总热量 $\int_0^1S_E\,dx=0.015$ 。与 Burgers 页的设计（源项的值与入口值）不同， $\mu_{\rm ref}$ 与 $Pr$ 出现在每个面的黏性通量里，所以对它们求导要经过所有的面（第 12 节）。页顶还可以改选“制造解”算例：入口相同，出口改为给定 $u=0.35$ 、 $T=T_{\rm out}$ （此时 $b_{\rm out}=T_{\rm out}$ ），并加上一个由解析解构造的固定源项；第 15 节用它测精度阶。两种算例是不同的边值问题（第 4 节）。

The design is the six-component vector $D=(a_0,a_1,\mu_{\rm ref},Pr,T_{\rm in},b_{\rm out})$ : two heating amplitudes, two transport parameters (the reference viscosity and the Prandtl number) and two boundary data, the inlet temperature $T_{\rm in}=D_4$ and the outlet datum $b_{\rm out}=D_5$ (here $p_b$ ). As is usual in the adjoint-optimization literature, the design variables are written $D$ . Section 5 makes $D$ the second argument of the discrete residual. The defaults are $D=(0.03,\,0,\,0.15,\,0.72,\,1,\,1)$ (bottom of Figure 1), so the heat added is $\int_0^1S_E\,dx=0.015$ . Unlike the design of the Burgers page (its source values and inlet value), $\mu_{\rm ref}$ and $Pr$ appear in the viscous flux of every face, so their derivatives run through all faces (Section 12). The top of the page also offers a manufactured case: the same inlet, an outlet that prescribes $u=0.35$ and $T=T_{\rm out}$ (then $b_{\rm out}=T_{\rm out}$ ), and a fixed source built from an analytic solution; Section 15 uses it to measure the order of accuracy. The two cases are different boundary-value problems (Section 4).

基准解是均匀状态附近的一个光滑小扰动（图 1）。没有热源时，均匀状态 $\rho=1$ 、 $u=0.35$ 、 $T=1$ 、 $p=1$ 满足方程和全部 5 个边界条件（ $T_{\rm in}=p_b=1$ ）。加上热源后，温度从入口的 $1$ 升到出口的 $1.0065$ ，约升高 $0.65\%$ 。压力几乎不变（与 $1$ 相差不超过 $0.15\%$ ）；在近似等压下受热，气体膨胀，密度下降约 $0.64\%$ ，而质量通量 $\rho u=0.35$ 不变，速度相应上升约 $0.65\%$ （以上数值取自连续参考解）。

The baseline solution is a smooth, small departure from a uniform state (Figure 1). Without heating, the uniform state $\rho=1$ , $u=0.35$ , $T=1$ , $p=1$ satisfies the equations and all 5 boundary conditions ( $T_{\rm in}=p_b=1$ ). With the heating, the temperature rises from $1$ at the inlet to $1.0065$ at the outlet, by about $0.65\%$ . The pressure hardly changes (it stays within $0.15\%$ of $1$ ); heated at nearly constant pressure, the gas expands: its density falls by about $0.64\%$ and, since the mass flux $\rho u=0.35$ is fixed, the velocity rises by about $0.65\%$ (values of the continuous reference).

对 $x$ 积分一次，三个分量各给出一个守恒关系（本算例动量方程没有源项）：

Integrating once in $x$ gives one conservation statement per component (the momentum equation has no source in this case):

\rho u=m_*=0.35,\qquad\rho u^2+p-\tau=\Pi,\qquad u(\rho E+p)-u\tau+q=\mathcal H(x)=\mathcal H(0)+\int_0^xS_E\,dx.

质量通量处处相同，总动量通量 $\Pi$ 处处相同，总能量通量 $\mathcal H$ 从入口起每经过一段就增加这一段加入的热量；附录 A 正是用 $\Pi$ 和 $\mathcal H$ 把方程组化成常微分方程。图 3把能量通量拆成三部分：对流 $u(\rho E+p)$ 、黏性功 $-u\tau$ 和导热 $q$ 。加入的热量 $0.015$ 中，约 $0.0082$ （55%）使出口处的对流通量高于入口，由流体从出口带走；约 $0.0071$ （47%）由热传导逆着流动方向从入口离开区间：入口温度固定为 $T_{\rm in}$ ，下游被加热，所以入口处 $T_x>0$ 、 $q=-\kappa T_x<0$ ；黏性功一项只变化 $-0.0003$ （−2%）。没有导热的 Euler 方程没有 $q$ ：加入的热量只能全部由流体从出口带走。格点页的圆点略偏离曲线：入口节点的三行被约束替换，节点 0 的半个控制体不再满足平衡，内部面上的质量通量是 $0.34984$ 而不是 $0.35$ 。这是 $N=12$ 时的离散误差，网格加密一倍后质量通量为 $0.34994$ （被替换的行见第 4 节）。

The mass flux is the same everywhere, so is the total momentum flux $\Pi$ , and the total energy flux grows along the duct by the heat added on the way; Appendix A uses exactly $\Pi$ and $\mathcal H$ to turn the equations into ordinary differential equations. Figure 3 splits the energy flux into three parts: convection $u(\rho E+p)$ , viscous work $-u\tau$ and conduction $q$ . Of the heat added, $0.015$ , about $0.0082$ (55%) raises the convective flux at the outlet above that at the inlet and is carried out of the outlet by the flow; about $0.0071$ (47%) leaves the interval through the inlet by conduction, against the flow: the inlet temperature is held at $T_{\rm in}$ and the gas downstream is heated, so at the inlet $T_x>0$ and $q=-\kappa T_x<0$ . The viscous-work term changes by only $-0.0003$ (−2%). The Euler equations have no $q$ : there, all of the added heat would leave through the outlet with the flow. On this node page the dots sit slightly off the curves: the three rows of the inlet node are replaced by constraints, so the half volume of node 0 no longer balances, and the interior faces carry the mass flux $0.34984$ instead of $0.35$ . This is the discretisation error at $N=12$ ; on a mesh twice as fine the mass flux is $0.34994$ (Section 4 describes the replaced rows).

图 3 · 能量通量沿流道的分解（默认参数；纵轴减去入口的对流通量 1.246437）。紫色：总能量通量 u(ρE + p) − uτ + q，从入口到出口的增量正是加入的热量 ΔH = ∫₀¹ S_E dx = 0.015；蓝色：对流 u(ρE + p)；橙色：导热 q = −κT_x；灰色虚线：黏性功 −uτ。曲线来自附录 A 的连续参考解，圆点是本页离散解在各面上的数值通量（N = 12；格点页内部面的质量通量是 0.34984 而不是 0.35，圆点因此略偏离曲线）。入口处 q < 0：约 47% 的热量由导热逆着流动从入口离开，约 55% 由流体从出口带走，黏性功约占 −2%。

Figure 3 · The energy flux along the duct, split into its parts (default parameters; the vertical axis subtracts the inlet convective flux 1.246437). Violet: the total energy flux u(ρE + p) − uτ + q, whose rise from the inlet to the outlet is exactly the heat added, ΔH = ∫₀¹ S_E dx = 0.015; blue: convection u(ρE + p); orange: conduction q = −κT_x; grey dashes: viscous work −uτ. The curves come from the continuous reference of Appendix A, the dots are this page's numerical fluxes on its faces (N = 12; the interior faces of this node page carry the mass flux 0.34984 instead of 0.35, so the dots sit slightly off the curves). At the inlet q < 0: about 47% of the heat leaves through the inlet by conduction, against the flow, about 55% is carried out of the outlet by the flow, and the viscous work accounts for about −2%.

第 7 节给出设计的目的。默认目标是温度跟踪：让温度 $T$ 接近目标温度 $T^*(x)=1+0.04\sin(\pi x)$ ，这是一个反设计问题。另两个可选目标是出口温度和入口热流 $q(0)$ ；后者正是上面由导热从入口送出的那部分热流（为负值）。本页计算所选目标对全部 6 个设计分量的导数。

Section 7 gives the purpose of the design. The default objective is temperature tracking: bring the temperature $T$ close to the target $T^*(x)=1+0.04\sin(\pi x)$ , an inverse-design problem. The two other objectives are the outlet temperature and the inlet heat flux $q(0)$ , which is exactly the conducted heat flux leaving through the inlet above (a negative number). The page computes the derivatives of the chosen objective with respect to all 6 design components.

这个方程组的作用

What these equations are for

在流体力学中，Navier–Stokes 方程组在 Euler 方程上加入分子黏性和热传导，是空气这类牛顿流体中黏性、导热流动的基本模型。黏性让相邻的流体层交换动量，导热让热量沿温度降低的方向扩散。雷诺数 $Re$ 比较对流与黏性扩散，Péclet 数 $Pe=Re\cdot Pr$ 比较对流与导热；二者都很大时，黏性和导热只在薄层里起主要作用（边界层、激波内部），一维方程组的一个经典应用就是计算正激波的内部结构。本页的 $Re$ 、 $Pe$ 都是 1 的量级，黏性和导热遍及整个区间，流动保持光滑、亚声速。

In fluid dynamics the Navier–Stokes equations add molecular viscosity and heat conduction to the Euler equations; they are the basic model of viscous, heat-conducting flow of a Newtonian fluid such as air. Viscosity exchanges momentum between neighbouring layers of fluid, and conduction carries heat towards lower temperature. The Reynolds number $Re$ compares convection with viscous diffusion and the Péclet number $Pe=Re\cdot Pr$ compares it with conduction; when both are large, viscosity and conduction matter mainly in thin layers (boundary layers, the inside of a shock), and a classical use of the one-dimensional equations is the internal structure of a normal shock. On this page $Re$ and $Pe$ are of order one, so viscosity and conduction act over the whole interval and the flow stays smooth and subsonic.

在本课程里，本页排在拟一维 Euler 之后、二维层流 N–S 之前（图 4）。与前两页相比，它新增了下面几件事：

In this course the page comes after quasi-1-D Euler and before 2-D laminar N–S (Figure 4). Compared with the pages before it, it adds:

同一组未知量，多一个通量。每点仍是 $\rho$ 、 $\rho u$ 、 $\rho E$ 三个守恒量，与拟一维 Euler 相同；面上的通量从 $F^c$ 变成 $F^c-F^v$ 。
The same unknowns, one more flux. Each point still carries the three conservative quantities $\rho$ , $\rho u$ , $\rho E$ , as on the quasi-1-D Euler page; the face flux becomes $F^c-F^v$ instead of $F^c$ .
面上要用梯度。黏性通量需要 $u_x$ 、 $T_x$ 。离散时它们来自面两侧状态之差，而面物性来自两侧的平均状态：一个面有两条依赖路径（第 3 节），前向与伴随都要沿两条路径求导。
Gradients at faces. The viscous flux needs $u_x$ and $T_x$ . Discretely they come from differences of the states on the two sides of a face, while the face properties come from the averaged state: a face has two dependency paths (Section 3), and the tangent and the adjoint must follow both.
随温度变化的物性。选 Sutherland 物性时， $\mu(T)$ 与 $\kappa(T)$ 把对温度的导数带进每个黏性通量。
Temperature-dependent properties. With Sutherland transport, $\mu(T)$ and $\kappa(T)$ bring derivatives with respect to temperature into every viscous flux.
设计变量进入物理模型。 $\mu_{\rm ref}$ 与 $Pr$ 出现在每个面的黏性通量里，而不只在源项或边界上（第 12 节）。
Design variables inside the physics. $\mu_{\rm ref}$ and $Pr$ appear in the viscous flux of every face, not only in a source or at a boundary (Section 12).
更多的边界条件。入口 3 个、出口 2 个，出口条件本身含应力与热流；格点页用约束替换 5 行残差（第 4 节）。
More boundary conditions. 3 at the inlet and 2 at the outlet, where the conditions involve the stress and the heat flux; the node page replaces 5 residual rows by constraints (Section 4).
独立的参考解。这里没有 Burgers 页那样的闭式解；附录 A 用射击法解一个独立的连续边值问题，给出参考解和连续伴随，第 15 节用制造解测精度阶。
An independent reference. There is no closed form as on the Burgers page; Appendix A solves an independent continuous boundary-value problem by shooting for a reference solution and the continuous adjoint, and Section 15 measures the order of accuracy with a manufactured solution.

图 4把本页放在整个网站的路线上。上排是一维 Burgers 页画出的主线，每往右一步主要增加一类新内容，每点的未知量依次为 1、3、4、4、5 个：一维 Burgers、拟一维 Euler、二维 Euler、二维层流 N–S、二维 RANS–SA。本页在主线下方：从拟一维 Euler 出发，先在一维里加入黏性和导热，再到二维层流 N–S 加上第二个空间方向；主线上的二维 Euler 则先加方向、后加黏性。二维层流 N–S 页使用同样的物性模型（常物性或 Sutherland 黏度， $\kappa=c_p\mu/Pr$ ），但应力是 $2\times2$ 张量，面上的梯度要由最小二乘重构得到。图中每个方框都链接到该页的第 1 节（同一种格式）。

Figure 4 places this page on the route through the site. The top row is the main line drawn on the 1-D Burgers page; each step to the right mainly adds one kind of ingredient, and the unknowns per point are 1, 3, 4, 4, 5: 1-D Burgers, quasi-1-D Euler, 2-D Euler, 2-D laminar N–S, 2-D RANS–SA. This page sits below the main line: starting from quasi-1-D Euler it adds viscosity and conduction while staying in one dimension, and 2-D laminar N–S then adds the second direction; along the main line, 2-D Euler adds the direction first and viscosity afterwards. The 2-D laminar N–S page uses the same transport model (constant or Sutherland viscosity, $\kappa=c_p\mu/Pr$ ), but its stress is a $2\times2$ tensor and its face gradients come from least-squares reconstructions. Every box of the figure links to Section 1 of its page (same scheme).

图 4 · 本页在网站路线上的位置。上排是一维 Burgers 页的路线：每往右一步主要增加一类新内容，每点的未知量依次是 1、3、4、4、5 个。本页在下方：从拟一维 Euler 出发，在一维里加入黏性应力和热传导（每点 3 个未知量），再通向二维层流 N–S。未知量取自各页的默认设置；点击方框跳到该页的第 1 节（同一种格式）。

Figure 4 · Where this page sits on the route through the site. The top row is the route of the 1-D Burgers page: each step to the right mainly adds one kind of ingredient, and the unknowns per point are 1, 3, 4, 4, 5. This page sits below: it starts from quasi-1-D Euler, adds viscous stress and heat conduction while staying in one dimension (3 unknowns per point), and leads on to 2-D laminar N–S. The unknowns come from each page's default settings; click a box to go to Section 1 of that page (same scheme).

本页要回答的问题

The question this page answers

本页只回答一个问题：设计 $D$ 的 6 个分量各改变一点，目标 $J$ 改变多少？页首的全页路线图把后面各节串成一条线。按页首默认设置（二阶、 $N=12$ 、块 GS 预处理的 GMRES）实际计数，6 次前向求解共用 206 次 GMRES 迭代，得到 6 个方向导数；1 次伴随求解只用 35 次，就给出全部 6 个分量。

This page answers one question: how much does the objective $J$ change when each of the 6 components of the design $D$ changes a little? The page roadmap at the top strings the later sections together. Counted with the default settings at the top of the page (second order, $N=12$ , block-GS-preconditioned GMRES), the 6 forward solves take 206 GMRES iterations in total for 6 directional derivatives, while the single adjoint solve takes 35 and gives all 6 components.

要在计算机上回答这个问题，先要把连续方程组变成有限个未知量的代数方程。第 2 节选定网格，决定三个守恒量存放在哪里；第 3 节再把第 1 节的通量写到每个面上。

To answer it on a computer, the continuous equations must first become algebraic equations in finitely many unknowns. Section 2 chooses the mesh and decides where the three conservative quantities are stored; Section 3 writes the fluxes of Section 1 on every face.

第 2 节Section 2

网格：未知量与数组编号

The mesh: unknowns and array indexing

第 1 节的方程组在区间上的每一个 $x$ 都成立，计算机却只能存有限个数。本节选定网格，决定三个守恒量存放在哪里、数组怎样编号；第 3 节再把第 1 节的通量写到每个面上。

The equations of Section 1 hold at every $x$ of the interval, but a computer stores finitely many numbers. This section chooses the mesh, where the three conservative quantities are stored and how the arrays are indexed; Section 3 then writes the fluxes of Section 1 on every face.

网格、数组和索引

Mesh, arrays and indices

把 $[0,1]$ 分为 $N$ 个等长区间， $h=1/N$ 。存储点为 $x_i=ih$ ， $i=0,\ldots,N$ ，共有 $M=N+1$ 个状态块。内部对偶面是相邻节点中点，两端面为 0、1；端点控制体宽度为 $h/2$ ，内部为 $h$ 。端点状态也是未知量，边界条件通过残差行替换施加。

Split $[0,1]$ into $N$ equal intervals with $h=1/N$ . Store $M=N+1$ state blocks at $x_i=ih$ , $i=0,\ldots,N$ . Interior dual faces lie halfway between nodes and endpoint faces are 0 and 1. Endpoint volumes have width $h/2$ , interior volumes $h$ . Endpoint states remain unknowns; boundary conditions replace residual rows.

U_i=\begin{bmatrix}\rho_i\\m_i\\e_i\end{bmatrix},\quad m_i=(\rho u)_i,\quad e_i=(\rho E)_i,\quad \mathbf U_h=(U_0^T,\ldots,U_{M-1}^T)^T\in\mathbb R^{3M}.

代码从 0 开始编号， $U_i$ 的三个分量在全局数组中分别占 $3i,3i+1,3i+2$ ；公式里的分量下标写作 1、2、3（质量、动量、能量），对应槽位 0、1、2。这里 $m_i$ 是动量密度，不是网格点数； $e_i$ 是总能量密度，不是比内能。内部面 $f$ 的左、右状态编号是 $L=f-1$ 、 $R=f$ 。下图横线上的点是状态存储位置，竖线是控制体面；数值来自当前求解。

Code indices start at zero: the components of $U_i$ occupy global slots $3i,3i+1,3i+2$ ; in formulas the component subscripts 1, 2, 3 (mass, momentum, energy) correspond to slots 0, 1, 2. Here $m_i$ is momentum density, not a point count, and $e_i$ is total energy density, not specific internal energy. Interior face $f$ has left and right indices $L=f-1$ , $R=f$ . Dots below locate stored states and vertical lines mark control-volume faces; values come from the current solve.

u_i=\frac{m_i}{\rho_i},\quad E_i=\frac{e_i}{\rho_i},\quad T_i=\frac{E_i-u_i^2/2}{c_v},\quad p_i=\rho_iT_i,\quad c_i=\sqrt{\gamma T_i}.

格点页的存储值是节点守恒量，由它转换的温度是节点处的点值；精度检验与解析节点值比较（第 15 节）。下标 $h$ 标记整个离散问题，不表示连续场（例如连续伴随）的采样，也不是差分验证的参数步长。

On the node page the stored values are nodal conservative states, and the temperature converted from them is a point value at the node; accuracy checks compare with exact nodal values (Section 15). The subscript $h$ labels the discrete problem, not samples of a continuous field (such as a sampled continuous adjoint) or a parameter-difference step.

图 5 · 同一条计算链，前向传播扰动，反向累加伴随。Figure 5 · One computational chain: forward perturbations and reverse adjoint accumulation.

网格和数组定下之后，下一步是面上的通量：第 3 节由两侧存储的状态算出每个内部面上的数值通量，对流部分来自重构的左右状态和数值耗散，黏性部分来自相邻状态之差。图 5 预览了整条计算链：第 3 节起逐步搭建，第 8、9 节对它求导。

With the mesh and the arrays fixed, the next step is the flux on a face: Section 3 computes the numerical flux on every interior face from the stored states on its two sides, the convective part from reconstructed left and right states with numerical dissipation, the viscous part from differences of neighbouring states. Figure 5 previews the whole chain: Section 3 onwards builds it, and Sections 8 and 9 differentiate it.

第 3 节Section 3

面上的通量：重构与黏性梯度

The flux on a face: reconstruction and viscous gradients

本节实验：先读下表每个面的离散梯度、面物性、应力与热流，再对照图 6 的两条路径：哪几列来自面平均状态，哪几列来自左右原始量之差？

Experiment: Read the face gradients, face transport, stress and heat flux in the table below, then match them with the two paths of Figure 6: which columns come from the averaged face state and which from differences of left and right primitive values?

第 2 节决定了状态存放在哪里、数组怎样编号；本节从这些存储的状态出发，算出内部面上的通量：先分清物理模型与离散运算，再写对流面状态的重构和黏性面的两条插值路径。

Section 2 decided where the states are stored and how the arrays are indexed; this section starts from those stored states and computes the flux on an interior face: first the line between the physical model and the discrete operations, then the reconstruction of the convective face states and the two interpolation paths of a viscous face.

物理模型与离散运算的分界

Separate the physical model from discrete operations

面上用的物理通量就是第 1 节的对流通量与黏性通量。下面先给出状态函数；它们接收一个三分量状态，而不是整个网格。第 1 节本构关系中的速度和温度空间导数，在程序里必须由明确的差分模板替代；从本节起，黏性面梯度统一记作 $G_f^u,G_f^T$ 。

The physical flux on a face is the convective and viscous flux of Section 1. The state functions below accept one three-component state, not the whole mesh. Spatial derivatives in the continuous constitutive law of Section 1 must be replaced by specified difference stencils in the program. Discrete viscous face gradients are denoted $G_f^u,G_f^T$ throughout the main text.

F^c(U)=\begin{bmatrix}m\\mu+p\\u(e+p)\end{bmatrix},\quad F_f^v=\begin{bmatrix}0\\\tau_f\\u_f^v\tau_f-q_f\end{bmatrix},\quad \tau_f=\frac43\mu_f G_f^u,\quad q_f=-\kappa_fG_f^T.

这就是第 1 节的 $\tau=\frac43\mu u_x$ 与 $q=-\kappa T_x$ ，只是面上的 $u_x$ 、 $T_x$ 换成了差分 $G_f^u$ 、 $G_f^T$ 。总通量为对流减黏性：黏性能量通量是 $u\tau-q$ ，减去它后热流以 $+q_f$ 进入净通量；而 $q_f=-\kappa_fG_f^T$ 本身带负号，所以净通量中的导热项是 $-\kappa_fG_f^T$ 。

These are the $\tau=\frac43\mu u_x$ and $q=-\kappa T_x$ of Section 1 with $u_x$ and $T_x$ replaced on the face by the differences $G_f^u$ and $G_f^T$ . Net flux is convection minus viscosity: the viscous energy flux is $u\tau-q$ , so subtracting it puts $+q_f$ into the net flux, and since $q_f=-\kappa_fG_f^T$ itself carries a minus sign, the conduction term in the net flux is $-\kappa_fG_f^T$ .

\mu_f=\begin{cases}\mu_{\rm ref}&\text{constant},\\\mu_{\rm ref}\dfrac{1.4(T_f^v)^{3/2}}{T_f^v+0.4}&\text{Sutherland},\end{cases}\qquad \kappa_f=\frac{c_p}{Pr}\mu_f.

面物性按第 1 节的物性模型，在面温度 $T_f^v$ （下文）上求值。一阶／二阶开关只改变对流和目标求积状态的重构；两种模式都保留中心黏性差分、同一边界闭合和物性模型。格点端点的三点单侧梯度对光滑点值为二阶；格心页的边界模板作用于单元平均，不能据此宣称也为二阶（见格心页第 4 节）。Rusanov 形式的数值耗散不等于分子黏性。

The face transport follows the transport law of Section 1, evaluated at the face temperature $T_f^v$ (below). The order switch changes reconstruction for convection and objective quadrature states. Both modes retain centred viscous differences, the same boundary closures and the same transport law. Three-point node endpoint gradients are second order for smooth point values; the boundary stencils of the cell page act on cell averages and are not second order for that reason (Section 4 of the cell page). Rusanov-form numerical dissipation is not molecular viscosity.

把二阶重构写成有限个权重

Write second-order reconstruction as finite weights

U_i^{\rm rec}(x)=\sum_jW_{ij}(x)U_j,\qquad d_i(x)=\frac{x-x_i}{h}.

U_i^{\rm rec}(x)=\begin{cases}U_i&\text{order 1},\\U_i+\dfrac{d_i}2(U_{i+1}-U_{i-1})&\text{order 2, interior},\\(1-d_0)U_0+d_0U_1& i=0,\ \text{order 2},\\(1+d_{M-1})U_{M-1}-d_{M-1}U_{M-2}&i=M-1,\ \text{order 2}.\end{cases}

重构直接作用于守恒量，没有 limiter。 $W_{ij}$ 是当前固定网格上的已知数，不是待求的连续微分算子。在内部面取 $U_f^-=U_L^{\rm rec}(x_f)$ 、 $U_f^+=U_R^{\rm rec}(x_f)$ 。边界面另按第 4 节处理，不能把内部双侧公式套到不存在的外侧单元。

Reconstruction acts on conservative variables without a limiter. Each $W_{ij}$ is a known number on the fixed mesh, not an unknown continuous differential operator. At an interior face use $U_f^-=U_L^{\rm rec}(x_f)$ and $U_f^+=U_R^{\rm rec}(x_f)$ . Section 4 supplies boundary closures; do not apply an interior two-sided formula to a nonexistent exterior cell.

\widehat F_f^c=\frac{F^c(U_f^-)+F^c(U_f^+)}2-\frac{\alpha_f}2(U_f^+-U_f^-),\quad \alpha_f=\frac{u_f^-+c_f^-+u_f^++c_f^+}{2},\quad H_f=\widehat F_f^c-F_f^v.

当前开放参数限定在正向、光滑、亚声速分支。平均波速用于该分支的耗散，不宣称是任意间断下的最大特征速度界。即使一阶重构， $\alpha_f$ 仍依赖状态，精确求导必须保留其变化。

The exposed parameter range targets a smooth positive-velocity subsonic branch. This averaged speed supplies dissipation on that branch; it is not claimed to bound all characteristic speeds across arbitrary jumps. Even with first-order reconstruction, $\alpha_f$ remains state dependent and must be differentiated.

黏性面有两条不同的插值路径

Two distinct paths to a viscous face

U_f^v=\frac{U_L+U_R}{2},\quad u_f^v=u(U_f^v),\quad T_f^v=T(U_f^v),\qquad G_f^u=\frac{u_R-u_L}{h},\quad G_f^T=\frac{T_R-T_L}{h}.

面物性使用平均守恒状态转换出的温度；梯度使用左右原始变量之差。一般有 $T((U_L+U_R)/2)\ne(T(U_L)+T(U_R))/2$ 。这是两个不同的计算节点，tangent 与反向传播都必须保留两条路径。代码变量 ux、tx 对应这里的离散梯度，不代表解析求导。

Face transport uses temperature converted from the averaged conservative state, whereas gradients difference the separately converted primitive states. In general $T((U_L+U_R)/2)\ne(T(U_L)+T(U_R))/2$ . These are distinct computational nodes; both tangent and reverse propagation retain both paths. Code variables ux and tx store these discrete gradients, not analytic derivatives.

图 6 · 黏性面的两条离散依赖路径。上路先平均守恒量；下路先转换再作差。前向沿连线传播，反向把两路贡献累加回同一存储状态。Figure 6 · Two discrete viscous-face paths: average conservative states on the upper path, convert then difference on the lower path. Forward propagation follows the connections; reverse accumulation combines both paths at the stored states.

以上公式处理的是内部对偶面：两侧各有一个节点。端点节点本身也是未知量；两个端面用端点状态的物理通量，黏性梯度改用单侧差分；面循环结束后，入口 3 行和出口 2 行残差被边界约束替换。第 4 节给出这些单侧梯度和约束。

These formulas treat interior dual faces, with a node on each side. The endpoint nodes are unknowns too: the two end faces use the physical flux of the endpoint state and one-sided differences for the viscous gradients, and after the face loop 3 inlet rows and 2 outlet rows of the residual are replaced by boundary constraints. Section 4 gives these one-sided gradients and constraints.

第 4 节Section 4

边界条件：规定哪些量，怎样进入通量

Boundary conditions: what is prescribed and how it enters the fluxes

先列出物理约束，再写成实际边界映射

From physical constraints to implemented boundary maps

第 1 节从方程的阶与特征的方向数出：入口要 3 个条件，出口要 2 个。本节把这些条件写成实际的边界映射。

Section 1 counted, from the order of the equations and the direction of the characteristics, 3 conditions at the inlet and 2 at the outlet. This section writes them as the implemented boundary maps.

加热问题入口规定 $\rho_b=1,u_b=0.35,T_b=D_4$ ，出口规定法向总应力 $p_f-\tau_f=D_5$ 与 $q_f=0$ 。它没有同时指定出口压力和零应力。制造解问题仍使用同一入口，出口改为规定 $u_b=0.35,T_b=D_5$ 。两种闭合是不同的边值问题，不能把它们的伴随边界项混用。

The heated case prescribes $\rho_b=1,u_b=0.35,T_b=D_4$ at the inlet, and total normal stress $p_f-\tau_f=D_5$ with $q_f=0$ at the outlet. It does not simultaneously prescribe pressure and zero stress. The manufactured case uses the same inlet but prescribes outlet $u_b=0.35,T_b=D_5$ . These are different boundary-value problems and their adjoint boundary contributions cannot be interchanged.

格点：先计算通量，最后替换五行

Nodes: evaluate fluxes, then replace five rows

G_0^\theta=\frac{-3\theta_0+4\theta_1-\theta_2}{2h},\quad G_M^\theta=\frac{3\theta_{M-1}-4\theta_{M-2}+\theta_{M-3}}{2h},\qquad \theta\in\{u,T\}.

边界对流通量取端点存储状态的物理通量，黏性通量使用上面的单侧差分。随后，入口三个分量替换为以下约束（约束向量写作 $\zeta$ ；字母 $c$ 在本页只表示声速）。注意：存储的是守恒量，而约束包含速度、温度，因此这些约束的 Jacobian 不是三阶单位阵。

Boundary convection uses the physical flux of the endpoint stored state; viscosity uses the one-sided differences above. The inlet's three rows are then replaced by the following constraints (the constraint vector is written $\zeta$ ; on this page the letter $c$ means only the sound speed). Stored variables are conservative, while constraints include velocity and temperature, so their Jacobian is not a three-by-three identity.

\zeta_{\rm in}=\begin{bmatrix}\rho_0-1\\u(U_0)-0.35\\T(U_0)-D_4\end{bmatrix},\quad\zeta_{\rm out}=\begin{cases}(p_{M-1}-\tau_M-D_5,\ q_M)^T&\text{heat},\\(u(U_{M-1})-0.35,\ T(U_{M-1})-D_5)^T&\text{MMS}.\end{cases}

出口只替换动量行和能量行，质量平衡行保留。压力／应力约束依赖端点及单侧梯度涉及的两个上游点，热流约束同样具有邻点耦合。反向传播先提取这五行的种子，再把旧守恒行种子清零；约束反传可继续写入邻点，不能把所有边界状态导数清零。

Only outlet momentum and energy rows are replaced; the mass balance remains. The pressure/stress constraint couples the endpoint to the two upstream points in the one-sided gradient, as does the heat-flux constraint. Reverse propagation first extracts the five constraint seeds, then zeros the overwritten balance seeds. Constraint derivatives can still accumulate into neighbours; do not zero all boundary state derivatives.

下面先比较端点状态及应力／热流，再核对净边界通量减体源项与原始平衡行之和。格点版的最终约束行之和不是守恒收支。条件数量本身不能证明连续问题适定或离散矩阵非奇异；这些属于不同层次，需结合当前分支与验证结果判断。

The panel compares endpoint states, stress and heat flux, then checks boundary net flux minus volume sources against the sum of raw balance rows. For nodes, summing final constraint rows is not a conservation budget. Counting conditions alone proves neither continuous well-posedness nor discrete nonsingularity; these are separate questions about the current branch and its verification.

第 5 节Section 5

残差：面循环、源项求积与边界覆盖

Residual: face loop, source quadrature and boundary rows

实际求导对象是一个数组到数组的函数

The derivative acts on an array-to-array function

R_h:\mathbb R^{3M}\times\mathbb R^6\longrightarrow\mathbb R^{3M},\qquad D=(a_0,a_1,\mu_{\rm ref},Pr,T_{\rm in},b_{\rm out})^T.

其中 $a_0,a_1$ 是两个热源幅度， $T_{\rm in}=D_4$ 是入口温度， $b_{\rm out}=D_5$ 是出口边界数据：加热算例是出口法向总应力的给定值 $p_b$ （ $p_f-\tau_f=D_5$ ），制造解算例是出口温度 $T_{\rm out}$ （第 4 节）。页顶参数选单在两种算例中都称它为“出口参数”；第 12–14 节与附录 A 的实验表格按算例标为 $p_b$ 或 $T_{\rm out}$ ，附录的公式写 $D_5$ ，它们是同一个分量。

Here $a_0,a_1$ are the two heat-source amplitudes, $T_{\rm in}=D_4$ is the inlet temperature and $b_{\rm out}=D_5$ is the outlet boundary datum: the prescribed normal total stress $p_b$ in the heated case ( $p_f-\tau_f=D_5$ ), and the outlet temperature $T_{\rm out}$ in the manufactured case (Section 4). The parameter menu at the top calls it "Outlet parameter" in both cases, the lab tables of Sections 12–14 and Appendix A label it $p_b$ or $T_{\rm out}$ according to the case, and the appendix formulas write $D_5$ : all name the same component.

固定网格、阶数、物性模式和边界类型后，程序接收 $\mathbf U_h,D$ ，输出 $R_h(\mathbf U_h,D)$ 。本文的离散伴随针对这个函数求导，不把连续偏微分算子及其边界条件直接搬来当作矩阵转置。制造源项在建立当前模型时固定；改变参数作局部求导时不重新制造。

With mesh, order, transport mode and boundary type fixed, the program maps $\mathbf U_h,D$ to $R_h(\mathbf U_h,D)$ . The discrete adjoint differentiates this function; it does not substitute a continuous differential operator and its boundary conditions for a matrix transpose. Manufactured forcing is frozen when the current model is created, not regenerated during local parameter differentiation.

R_{i}^{\rm raw}=H_{i+1}-H_i-S_i^h,\quad S_i^h=\sum_{\pm}\omega_{i\pm}S(x_{i\pm},D),\quad x_{i\pm}=\frac{x_i^-+x_i^+}{2}\pm\frac{V_i}{2\sqrt3},\quad\omega_{i\pm}=\frac{V_i}{2}.

$x_i^-,x_i^+$ 是控制体左右面， $V_i=x_i^+-x_i^-$ 。源项使用两点 Gauss 求积；这里写的是程序实际求和，不是未经离散的积分。残差没有除以体积。一个内部面只求值一次，向左加、向右减；这样成对消去来自装配本身，而不是两个独立通量近似相等。

The control-volume faces are $x_i^-,x_i^+$ and $V_i=x_i^+-x_i^-$ . Sources use two-point Gauss quadrature: this is the actual program sum, not an undiscretized integral. Residuals are not divided by volume. Compute each interior face once, adding to the left and subtracting from the right. Pairwise cancellation follows from assembly itself, not approximate equality of independently evaluated fluxes.

S(x,D)=S_{\rm fixed}(x)+e_3\left[a_0\sin^2(\pi x)+a_1\sin^2(2\pi x)\right],\qquad e_3=(0,0,1)^T.

加热模式的固定源项为零；制造解模式的固定源项由独立解析函数计算。它不依赖当前离散状态，所以本模型 $\partial S_i^h/\partial\mathbf U_h=0$ ；对热源参数的偏导则非零。将来若引入状态相关源项，这个结论不再成立。

The fixed source is zero in heating mode and comes from independent analytic functions in manufactured mode. It is independent of the current discrete state, so $\partial S_i^h/\partial\mathbf U_h=0$ in this model, while heat-parameter derivatives remain nonzero. This conclusion would change for a state-dependent source.

R_h=M_{\rm bc}R^{\rm raw}+E_{\rm bc}\zeta(\mathbf U_h,D).

$M_{\rm bc}$ 在五条被替换行取 0，其余对角元取 1； $E_{\rm bc}$ 把第 4 节的五个约束（向量 $\zeta$ ）插入相应行。这个写法把“覆盖数组元素”也纳入数学函数定义。

The diagonal mask $M_{\rm bc}$ is zero on the five replaced rows and one elsewhere; $E_{\rm bc}$ inserts Section 4's five constraints (the vector $\zeta$ ). Thus overwriting array entries is part of the mathematical function definition.

Rraw[:] = 0
for f = 0..M:
    traces = reconstruct_or_close_boundary(U, D, f)
    face_state, Gu, GT = viscous_face_data(U, D, f)
    H[f] = convective_flux(traces) - viscous_flux(face_state, Gu, GT, D)
    if f > 0: Rraw[f-1] += H[f]
    if f < M: Rraw[f]   -= H[f]
for i = 0..M-1: Rraw[i] -= two_point_source_quadrature(i, D)
R = copy(Rraw)
if node_centred: overwrite_five_boundary_rows(R, U, D)

播放器按这一装配顺序更新数组。先观察一个面的左右等量反号贡献，再跳到末尾检查源项及边界覆盖。原始量表、局部通量表和最终残差数组分别是不同层次，不能把局部通量的导数直接当作全局残差导数。

The player updates arrays in this assembly order. Inspect the equal and opposite contribution of one face, then advance to sources and boundary overwrites. Primitive states, local face fluxes and the final residual array are different levels; a local flux derivative alone is not the global residual derivative.

第 6 节Section 6

Jacobian 与块 GS：精确导数与近似求解矩阵

Jacobian and block GS: the exact derivative and the approximate solver matrix

精确 Jacobian 与近似求解矩阵

Exact Jacobian and approximate solve matrix

A_h=\frac{\partial R_h}{\partial\mathbf U_h}\in\mathbb R^{3M\times3M},\qquad C_h=\frac{\partial R_h}{\partial D}\in\mathbb R^{3M\times6}.

$A_h$ 包含重构、原始变量转换、数值耗散、黏性差分、温度物性和边界闭合的全部状态依赖。 $C_h$ 则固定状态数组，只改变参数。目标函数不属于残差 Jacobian；它另行给出第 7 节的右端向量。二阶重构使单个面的导数访问更多邻点，因此即使一维网格也不能想当然地只保留最近邻块。

The exact $A_h$ includes every state dependency through reconstruction, primitive conversion, numerical dissipation, viscous differences, temperature-dependent transport and boundary closure. The parameter Jacobian $C_h$ changes parameters while holding the state array fixed. The objective is not part of the residual Jacobian; it supplies the separate right-hand side in Section 7. Second-order reconstruction extends face dependencies to more neighbours, so even in one dimension nearest-neighbour blocks alone are not generally exact.

P_h=\left.\frac{\partial R_h^{(1)}}{\partial\mathbf U_h}\right|_{\mathbf U_h},\quad P_h=\operatorname{diag}(P_h)+L_P+U_P,\qquad B_h=(\operatorname{diag}(P_h)+L_P)^{-1}.

代码把同一当前状态送入一阶对流版本，构造 $P_h$ ；黏性、物性和边界仍完整保留。 $\operatorname{diag}(P_h),L_P,U_P$ 分别是块对角、严格块下三角和严格块上三角，每个块为 $3\times3$ 。 $B_h$ 表示从零初值出发的一次前向块 GS 扫描，不是 $P_h^{-1}$ 。即使一阶模式 $P_h=A_h$ ，一次扫描通常也不是线性方程的精确解。

The code constructs $P_h$ from first-order convection at the same current state, retaining full viscosity, transport and boundaries. Its block diagonal, strictly lower and strictly upper parts are $\operatorname{diag}(P_h),L_P,U_P$ , with $3\times3$ blocks. The map $B_h$ is one forward block GS sweep starting from zero, not $P_h^{-1}$ . Even when first-order mode gives $P_h=A_h$ , one sweep is generally not the exact linear solution.

(B_hr)_i=(P_{ii})^{-1}\left(r_i-\sum_{j<i}P_{ij}(B_hr)_j\right),\quad i=0,\ldots,M-1.

先解非线性原始问题

Solve the nonlinear primal problem first

A_h(\mathbf U_h^{(k)})\Delta\mathbf U_h=-R_h(\mathbf U_h^{(k)},D),\qquad \mathbf U_h^{(k+1)}=\mathbf U_h^{(k)}+\beta_k\Delta\mathbf U_h.

这是 Newton 修正；每一步重新线性化，用右预处理 GMRES 解线性系统，然后回溯选择步长。实际代码要求重构面状态保持正密度、正温度及指定正向亚声速范围，同时残差满足充分下降。改变步长是找解的策略，不是额外物理项。稳态伴随在最终收敛状态重新建立精确 Jacobian，不对回溯次数求导。

This is a Newton correction. Each iteration relinearizes, solves with right-preconditioned GMRES, and backtracks the step. The implementation checks reconstructed face states for positive density and temperature and the prescribed positive subsonic range, together with sufficient residual decrease. Step control is a solution strategy, not an extra physical term. The steady adjoint rebuilds the exact Jacobian at the converged state; it does not differentiate the number of backtracking steps.

下面第一张图是 Newton 非线性残差；选择某次迭代可查看对应状态。矩阵图为教学检查额外显式组装，不能据此把主导数接口称为非 matrix-free：主接口仍通过计算图计算乘积，显式矩阵用于结构观察、预处理及附录 C 的直接法。

The first plot shows the nonlinear Newton residual; selecting an iteration displays its state. The matrix plot is an extra explicit assembly for inspection. Main derivative interfaces still evaluate graph products without storing the exact matrix; explicit matrices serve visualization, preconditioning and Appendix C's direct method.

第 7 节Section 7

目标函数：温度跟踪、出口温度与入口热流

Objectives: temperature tracking, outlet temperature and inlet heat flux

先离散目标，再形成偏导向量

Discretize the objective before forming derivative vectors

J_h:\mathbb R^{3M}\times\mathbb R^6\to\mathbb R,\quad j_U=\left(\frac{\partial J_h}{\partial U_{h,0}},\ldots,\frac{\partial J_h}{\partial U_{h,3M-1}}\right)^T,\quad j_D=\left(\frac{\partial J_h}{\partial D_0},\ldots,\frac{\partial J_h}{\partial D_5}\right)^T.

本页统一把 $j_U,j_D$ 定义为列向量，避免把偏导数行与梯度列直接相加。求 $j_U$ 固定参数；求 $j_D$ 固定所有状态槽位。伴随右端是 $-j_U$ ，而最终六维设计梯度还包含 $j_D$ 。它们的维数及物理含义不同。

This page consistently defines $j_U,j_D$ as columns, avoiding addition of derivative rows to gradient columns. Hold parameters fixed for $j_U$ , and every state slot fixed for $j_D$ . The adjoint right-hand side is $-j_U$ , whereas the final six-dimensional design gradient also includes $j_D$ . Their dimensions and meanings are different.

温度跟踪：从重构到求积

Temperature tracking: reconstruction through quadrature

\Theta_{i\pm}=T(U_i^{\rm rec}(x_{i\pm})),\quad \Theta^*_{i\pm}=1+0.04\sin(\pi x_{i\pm}),\qquad J_{T,h}=\frac12\sum_i\sum_{\pm}\omega_{i\pm}(\Theta_{i\pm}-\Theta^*_{i\pm})^2.

求积点与权重和源项一致，但求积状态需要重构，随后才转换成温度。固定网格和固定解析目标意味着点、权重与目标值对设计的 tangent 均为零；本页不做几何设计。用 $\dot U_j$ 表示存储状态的 tangent（前向模式的方向导数，点号约定见第 8 节），则逐项链式法则给出：

Quadrature points and weights match the source quadrature, but objective states are reconstructed before conversion to temperature. With this fixed mesh and fixed analytic target, points, weights and target samples have zero tangents with respect to the design; this page does not optimize geometry. With $\dot U_j$ the tangent of a stored state (a forward-mode directional derivative; Section 8 sets out the dot convention), the chain rule gives:

\dot J_{T,h}=\sum_{i,\pm}\omega_{i\pm}(\Theta_{i\pm}-\Theta^*_{i\pm})\,T_U(U_i^{\rm rec}(x_{i\pm}))\sum_jW_{ij}(x_{i\pm})\dot U_j.

(j_U)_j=\sum_{i,\pm}W_{ij}(x_{i\pm})\,\omega_{i\pm}(\Theta_{i\pm}-\Theta^*_{i\pm})\,T_U(U_i^{\rm rec}(x_{i\pm}))^T,\qquad j_D=0\quad\text{(tracking, fixed state)}.

这里 $(j_U)_j$ 表示三个分量组成的块。求积权重已经进入右端，伴随方程中不能再乘一次体积。二阶目标会把一个求积点的权重分发给多个存储状态；只在最近状态上写目标梯度会漏掉重构路径。

Here $(j_U)_j$ is a three-component block. Quadrature weights already enter the right-hand side and must not be multiplied in again. A second-order objective scatters one quadrature point's contribution to several stored states; writing it only to the nearest state omits the reconstruction path.

端点目标：显式参数项不能省

Endpoint objectives retain explicit parameter terms

J_{{\rm out},h}=T_f\big|_{f=M},\qquad J_{q,h}=q_0=-\kappa_0G_0^T,\qquad \dot J_{q,h}=-\kappa_0\dot G_0^T-G_0^T\dot\kappa_0.

热流目标即使固定状态，也会随 $\mu_{\rm ref},Pr$ 改变；格心入口还显式依赖规定温度。制造解模式的出口温度在最终收敛后由边界给定，对热源幅度的总导数应为零。格心版可以直接通过目标的参数项得到，格点版则经强边界约束得到同一结论；不能把两版的中间偏导项强行要求相等。

Even with fixed state, heat flux changes with $\mu_{\rm ref},Pr$ ; the cell inlet also explicitly depends on prescribed temperature. In manufactured mode the converged outlet temperature is prescribed, so its total heat-source derivative should vanish. Cells can express this directly through the objective's parameter term; nodes reach the same conclusion through their strong boundary constraint. Intermediate partial derivatives need not match between the two schemes.

第 8 节Section 8

前向：沿残差依赖链传播方向导数

Tangent: propagate directional derivatives along the residual chain

本节实验：一次 tangent 计算给出残差对指定方向的变化，不是已经求出了状态对设计的敏感度。比较第 10 节的线性求解。

Experiment: One tangent application gives the residual change along a supplied direction. It is not yet the solved state sensitivity; compare the linear solve in Section 10.

点号表示有限维方向导数

A dot denotes a finite-dimensional directional derivative

\dot R_h=\left.\frac{d}{d\epsilon}R_h(\mathbf U_h+\epsilon v,D+\epsilon\dot D)\right|_{\epsilon=0}=A_hv+C_h\dot D,\qquad \dot J_h=j_U^Tv+j_D^T\dot D.

状态种子 $v\equiv\dot{\mathbf U}_h$ 是存储状态的 tangent， $\dot D$ 是六维参数种子（设计方向）。本页按自动微分的记法书写离散程序的导数：点号是前向模式（tangent）导数，第 9 节的上横线是反向模式（adjoint）导数；变分符号 $\delta$ 只用于附录 A 中连续函数的变分。这里没有时间积分，点号也不是时间导数。计算 $A_hv$ 时令 $\dot D=0$ ；计算 $C_h\dot D$ 时令 $v=0$ 。下面沿实际残差依赖链逐步传播，而不是近似差分一个 Jacobian。

The state seed $v\equiv\dot{\mathbf U}_h$ is the tangent of the stored states and $\dot D$ is the six-dimensional parameter seed (the design direction). This page writes derivatives of the discrete program in the notation of algorithmic differentiation: a dot is a forward-mode (tangent) derivative and the overbar of Section 9 a reverse-mode (adjoint) derivative, while the variation symbol $\delta$ is used only for variations of continuous functions in Appendix A. There is no time integration, so a dot is not a time derivative either. Set $\dot D=0$ for $A_hv$ and $v=0$ for $C_h\dot D$ . Follow the actual residual dependency chain below rather than approximating a Jacobian by finite differences.

第一步：守恒量转换的局部 Jacobian

Step 1: the local conservative-to-primitive Jacobian

\dot U=(\dot\rho,\dot m,\dot e)^T,\quad \dot u=\frac{\dot m-u\dot\rho}{\rho},\quad\dot T=\frac{\dot e-E\dot\rho}{\rho c_v}-\frac{u\dot u}{c_v},\quad\dot p=T\dot\rho+\rho\dot T,\quad\dot c=\frac{\gamma}{2c}\dot T.

u_U=\begin{bmatrix}-u/\rho&1/\rho&0\end{bmatrix},\qquad T_U=\frac1{\rho c_v}\begin{bmatrix}u^2-E&-u&1\end{bmatrix}.

上式 $\dot U$ 是一个局部三分量状态的 tangent，顺序为密度、动量密度、能量密度；对应代码分量槽位 0、1、2。网格位置另用 $U_i$ 的下标 $i$ 表示。 $u_U,T_U$ 是 $1\times3$ 行，分别与三分量 tangent 相乘。所有系数在当前基态求值。

Here $\dot U$ is the tangent of one local three-component state, ordered as density, momentum density and energy density, matching code component slots 0, 1 and 2. The subscript $i$ in $U_i$ separately identifies a mesh location. The rows $u_U,T_U$ have size $1\times3$ and multiply a three-component tangent. Evaluate every coefficient at the current base state.

气体常数为 1， $c_v=1/(\gamma-1)$ ，所以压力可以直接写成守恒量的函数。它的行导数 $p_U$ 与上式的 $\dot p=T\dot\rho+\rho\dot T$ 相同，后面的通量 Jacobian 与格点出口约束都用它：

With the gas constant 1 and $c_v=1/(\gamma-1)$ , the pressure is a function of the conservative variables alone. Its row $p_U$ is the same as $\dot p=T\dot\rho+\rho\dot T$ above, and the flux Jacobian and the node outlet constraint below use it:

p=\rho T=(\gamma-1)\Bigl(e-\frac{m^2}{2\rho}\Bigr),\qquad p_U=(\gamma-1)\begin{bmatrix}\tfrac12u^2&-u&1\end{bmatrix},\qquad c_U=\frac{\gamma}{2c}T_U.

第二步：对流面重构与耗散

Step 2: convective traces and dissipation

\dot U_f^- =\sum_jW_{Lj}(x_f)v_j,\quad\dot U_f^+=\sum_jW_{Rj}(x_f)v_j,\qquad\dot\alpha_f=\frac{\dot u_f^-+\dot c_f^-+\dot u_f^++\dot c_f^+}{2}.

\dot F^c(U)=\begin{bmatrix}\dot m\\u\dot m+m\dot u+\dot p\\(e+p)\dot u+u(\dot e+\dot p)\end{bmatrix},

\dot{\widehat F}_f^c=\frac{\dot F^c(U_f^-)+\dot F^c(U_f^+)}2-\frac{\dot\alpha_f}2(U_f^+-U_f^-)-\frac{\alpha_f}2(\dot U_f^+-\dot U_f^-).

冻结耗散速度会删去倒数第二项 $-\frac{\dot\alpha_f}2(U_f^+-U_f^-)$ ，只保留耗散部分的最后一项，得到的是另一套线性算子。它可能通过自己的转置点积测试，但并不是当前残差的精确导数。边界面使用物理通量，没有内部 Rusanov 耗散；必须先求导边界状态映射。

Freezing dissipation speed drops the penultimate term and produces a different linear operator. That operator may pass its own transpose test while failing to differentiate the actual residual. Boundary faces use physical fluxes without interior Rusanov dissipation; differentiate their state closure first.

对流通量 Jacobian

The convective flux Jacobian

把 $\dot F^c(U)$ 的三行按 $\dot\rho,\dot m,\dot e$ 合并，就是 Euler 通量的 Jacobian 矩阵（ $h_{\rm t}$ 是比总焓，气体常数为 1 时 $p/\rho=T$ ；本页的 $H_f$ 是面上的净通量，与它无关）：

Collecting the three rows of $\dot F^c(U)$ by $\dot\rho,\dot m,\dot e$ gives the Jacobian of the Euler flux ( $h_{\rm t}$ is the specific total enthalpy; with the gas constant 1, $p/\rho=T$ ; the $H_f$ of this page is the net flux through a face and has nothing to do with it):

\dot F^c(U)=F^c_U\dot U,\qquad F^c_U=\begin{bmatrix}0&1&0\\-\frac{3-\gamma}{2}u^2&(3-\gamma)u&\gamma-1\\u\bigl(\frac{\gamma-1}{2}u^2-h_{\rm t}\bigr)&h_{\rm t}-(\gamma-1)u^2&\gamma u\end{bmatrix},\qquad h_{\rm t}=\frac{e+p}{\rho}=E+T.

第二行来自 $\rho u^2+p=\frac{3-\gamma}{2}\frac{m^2}{\rho}+(\gamma-1)e$ ，第三行来自 $u(e+p)=\gamma\frac{me}{\rho}-\frac{\gamma-1}{2}\frac{m^3}{\rho^2}$ ，逐项对 $\rho,m,e$ 求偏导即得。例如在入口规定状态 $\rho=1$ 、 $u=0.35$ 、 $T=1$ ：

The second row comes from $\rho u^2+p=\frac{3-\gamma}{2}\frac{m^2}{\rho}+(\gamma-1)e$ and the third from $u(e+p)=\gamma\frac{me}{\rho}-\frac{\gamma-1}{2}\frac{m^3}{\rho^2}$ , differentiated term by term with respect to $\rho,m,e$ . At the prescribed inlet state $\rho=1$ , $u=0.35$ , $T=1$ , for example:

F^c_U=\begin{bmatrix}0&1&0\\-0.098&0.56&0.4\\-1.238&3.512&0.49\end{bmatrix},\qquad\text{eigenvalues }u-c,\ u,\ u+c=-0.833,\ 0.350,\ 1.533.

这三个特征值就是第 1 节图 2 的入口特征速度。第 9 节的对流转置直接用 $F^c_U(U_f^\pm)^T$ ；Rusanov 耗散还要 $u_U$ 与 $c_U$ ，它们都在两侧重构状态上求值。

These eigenvalues are the inlet characteristic speeds of Figure 2 in Section 1. The convective transpose of Section 9 uses $F^c_U(U_f^\pm)^T$ directly; the Rusanov dissipation also needs $u_U$ and $c_U$ , both evaluated on the two reconstructed traces.

第三步：两条黏性路径和物性

Step 3: both viscous paths and transport

\dot U_f^v=\frac{v_L+v_R}{2},\quad\dot u_f^v=u_U(U_f^v)\dot U_f^v,\quad\dot T_f^v=T_U(U_f^v)\dot U_f^v,

\dot G_f^u=\frac{u_U(U_R)v_R-u_U(U_L)v_L}{h},\quad\dot G_f^T=\frac{T_U(U_R)v_R-T_U(U_L)v_L}{h}.

\dot\mu_f=\mu_{T,f}\dot T_f^v+\frac{\mu_f}{\mu_{\rm ref}}\dot\mu_{\rm ref},\quad\mu_{T,f}=\begin{cases}0&\text{constant},\\\mu_f\left(\dfrac{3}{2T_f^v}-\dfrac1{T_f^v+0.4}\right)&\text{Sutherland},\end{cases}

\dot\kappa_f=\frac{c_p}{Pr}\dot\mu_f-\frac{\kappa_f}{Pr}\dot{Pr},\quad\dot\tau_f=\frac43(\mu_f\dot G_f^u+G_f^u\dot\mu_f),\quad\dot q_f=-\kappa_f\dot G_f^T-G_f^T\dot\kappa_f.

\dot F_f^v=\begin{bmatrix}0\\\dot\tau_f\\\tau_f\dot u_f^v+u_f^v\dot\tau_f-\dot q_f\end{bmatrix},\qquad\dot H_f=\dot{\widehat F}_f^c-\dot F_f^v.

第四步：黏性通量是面状态与面梯度的函数

Step 4: the viscous flux as a function of the face state and the face gradient

第三步逐项写出了各个 tangent。按输入分组，本构黏性通量只依赖三类量：面平均状态 $U_f^v$ （经过 $u_f^v$ 、 $T_f^v$ 以及物性 $\mu(T)$ 、 $\kappa(T)$ ）、两个面梯度 $G_f^u,G_f^T$ ，以及设计分量 $\mu_{\rm ref},Pr$ 。所以它的 tangent 是四个分块 Jacobian 之和（ $u_U,T_U$ 在 $U_f^v$ 上求值）：

Step 3 wrote every tangent term by term. Grouped by input, the constitutive viscous flux depends on three kinds of quantity only: the averaged face state $U_f^v$ (through $u_f^v$ , $T_f^v$ and the transport $\mu(T)$ , $\kappa(T)$ ), the two face gradients $G_f^u,G_f^T$ , and the design components $\mu_{\rm ref},Pr$ . Its tangent is therefore a sum of four block Jacobians ( $u_U,T_U$ evaluated at $U_f^v$ ):

F_f^v=\begin{bmatrix}0\\\tfrac43\mu_fG_f^u\\\tfrac43\mu_fu_f^vG_f^u+\kappa_fG_f^T\end{bmatrix},\qquad\dot F_f^v=\frac{\partial F^v}{\partial U^v}\dot U_f^v+\frac{\partial F^v}{\partial G^u}\dot G_f^u+\frac{\partial F^v}{\partial G^T}\dot G_f^T+\frac{\partial F^v}{\partial\mu_{\rm ref}}\dot\mu_{\rm ref}+\frac{\partial F^v}{\partial Pr}\dot{Pr},

\frac{\partial F^v}{\partial U^v}=\begin{bmatrix}0\\\tfrac43G_f^u\mu_{T,f}\,T_U\\\tau_f\,u_U+\bigl(\tfrac43u_f^vG_f^u\mu_{T,f}+G_f^T\kappa_{T,f}\bigr)T_U\end{bmatrix},\qquad\frac{\partial F^v}{\partial G^u}=\begin{bmatrix}0\\\tfrac43\mu_f\\\tfrac43\mu_fu_f^v\end{bmatrix},\qquad\frac{\partial F^v}{\partial G^T}=\begin{bmatrix}0\\0\\\kappa_f\end{bmatrix},

\frac{\partial F^v}{\partial\mu_{\rm ref}}=\frac{F_f^v}{\mu_{\rm ref}},\qquad\frac{\partial F^v}{\partial Pr}=\begin{bmatrix}0\\0\\q_f/Pr\end{bmatrix},\qquad\kappa_{T,f}=\frac{c_p}{Pr}\mu_{T,f}.

$\mu_f$ 与 $\kappa_f$ 都和 $\mu_{\rm ref}$ 成正比，所以整个本构黏性通量与 $\mu_{\rm ref}$ 成正比，对它的导数就是通量本身除以 $\mu_{\rm ref}$ ； $Pr$ 只出现在 $\kappa_f=c_p\mu_f/Pr$ 里，只改变能量分量中的导热。第 12 节用这两个恒等式把设计导数写成逐面求和。

Both $\mu_f$ and $\kappa_f$ are proportional to $\mu_{\rm ref}$ , so the whole constitutive viscous flux is proportional to $\mu_{\rm ref}$ and its derivative is the flux itself divided by $\mu_{\rm ref}$ ; $Pr$ appears only in $\kappa_f=c_p\mu_f/Pr$ and changes only the conduction of the energy component. Section 12 uses these two identities to write the design derivatives as sums over faces.

物性对温度的导数与热流

The temperature derivative of the transport, and the heat flux

Sutherland 律取对数后逐项求导：

Take the logarithm of the Sutherland law and differentiate term by term:

\ln\mu=\ln\mu_{\rm ref}+\ln\,1.4+\tfrac32\ln\,T-\ln(T+0.4)\quad\Longrightarrow\quad\frac{\mu_T}{\mu}=\frac{3}{2T}-\frac{1}{T+0.4},\qquad\frac{\mu_T}{\mu}\Big|_{T=1}=\frac32-\frac1{1.4}=\frac{11}{14}.

所以在 $T=1$ 附近温度升高 1% 时黏度约升高 0.79%；常物性时 $\mu_T=0$ 。热流 $q_f=-\kappa_fG_f^T$ 同时经过梯度与导热系数，导热系数又经过温度和两个设计分量：

Near $T=1$ a 1% rise in temperature thus raises the viscosity by about 0.79%; with constant transport $\mu_T=0$ . The heat flux $q_f=-\kappa_fG_f^T$ runs through both the gradient and the conductivity, and the conductivity through the temperature and two design components:

\dot q_f=-\kappa_f\dot G_f^T-G_f^T\Bigl(\kappa_{T,f}\dot T_f^v+\frac{\kappa_f}{\mu_{\rm ref}}\dot\mu_{\rm ref}-\frac{\kappa_f}{Pr}\dot{Pr}\Bigr).

面梯度是一个有限模板

A face gradient is a finite stencil

内部面与两个端面的梯度都是存储状态的线性组合：

The gradients of interior faces and of the two end faces are linear combinations of stored states:

G_f^\theta=\sum_jw_{fj}\,\theta(U_j),\qquad\dot G_f^\theta=\sum_jw_{fj}\,\theta_U(U_j)\,v_j,\qquad\theta\in\{u,T\},

\text{interior: }w_{f,L}=-\tfrac1h,\ w_{f,R}=\tfrac1h;\qquad\text{inlet: }(w_{00},w_{01},w_{02})=\bigl(-\tfrac3{2h},\tfrac2h,-\tfrac1{2h}\bigr);\qquad\text{outlet: }(w_{M,M-1},w_{M,M-2},w_{M,M-3})=\bigl(\tfrac3{2h},-\tfrac2h,\tfrac1{2h}\bigr).

一阶与二阶的黏性模板相同。模板的 tangent 是“收集”：每个模板点贡献 $w_{fj}\theta_U(U_j)v_j$ ；第 9 节会看到它的伴随是用同一组权重“分发”（图 9）。

The viscous stencils are the same at first and second order. The tangent of a stencil gathers: every stencil point contributes $w_{fj}\theta_U(U_j)v_j$ ; Section 9 shows that its adjoint scatters with the same weights (Figure 9).

第五步：边界面、边界行与源项

Step 5: boundary faces, boundary rows and sources

以上内部面梯度在边界换成第 4 节的单侧线性组合，最后再用五条约束的 tangent 覆盖相应残差行。例如入口温度行的 tangent 是 $T_U(U_0)v_0-\dot D_4$ ，不是能量槽位的 tangent 减 $\dot D_4$ 。

At boundaries replace the interior gradient formulas by Section 4's one-sided combinations, then overwrite the five rows with the tangents of the constraints. For example the tangent of the inlet temperature row is $T_U(U_0)v_0-\dot D_4$ , not the energy-slot tangent minus $\dot D_4$ .

格点页的两个端面在节点上：对流通量取端点状态的物理通量（没有 Rusanov 耗散），梯度用第四步的三点单侧模板，tangent 照第一到四步写出。随后五行被约束覆盖，它们的 tangent 也是覆盖而不是累加：

The two end faces of the node page sit on nodes: the convective flux is the physical flux of the end state (no Rusanov dissipation), the gradients use the three-point one-sided stencils of Step 4, and the tangents follow Steps 1–4. The five rows are then overwritten by the constraints, and so are their tangents, which replace rather than accumulate:

\dot\zeta_{\rm in}=\begin{bmatrix}v_{0,1}\\u_U(U_0)v_0\\T_U(U_0)v_0-\dot D_4\end{bmatrix},\qquad\dot\zeta_{\rm out}=\begin{cases}\bigl(p_U(U_{M-1})v_{M-1}-\dot\tau_M-\dot D_5,\ \dot q_M\bigr)^T&\text{heat},\\\bigl(u_U(U_{M-1})v_{M-1},\ T_U(U_{M-1})v_{M-1}-\dot D_5\bigr)^T&\text{MMS},\end{cases}

\dot\tau_M=\tfrac43\bigl(\mu_M\dot G_M^u+G_M^u\dot\mu_M\bigr),\qquad\dot q_M=-\kappa_M\dot G_M^T-G_M^T\dot\kappa_M,\qquad\dot R_{h,r}=\dot\zeta_r\quad\text{on the five rows }r.

这里 $v_{0,1}$ 是 $v_0$ 的密度分量， $\mu_M,\kappa_M$ 在节点 $M-1$ 的温度上求值， $\dot G_M^\theta$ 用出口模板。出口只覆盖动量行与能量行，质量行仍是平衡行，仍收到 $\dot H_M$ 的质量分量。

Here $v_{0,1}$ is the density component of $v_0$ , $\mu_M,\kappa_M$ are evaluated at the temperature of node $M-1$ , and $\dot G_M^\theta$ uses the outlet stencil. The outlet overwrites only the momentum and energy rows; its mass row stays a balance row and still receives the mass component of $\dot H_M$ .

源项不依赖状态，它的 tangent 只来自两个热源幅度（ $Q_{ik}$ 是第 12 节的求积权重）：

The source does not depend on the state; its tangent comes from the two heat amplitudes only ( $Q_{ik}$ are the quadrature weights of Section 12):

\dot S_i^h=e_3\bigl(Q_{i0}\dot D_0+Q_{i1}\dot D_1\bigr),\qquad\dot R^{\rm raw}_i\;{-}\!{=}\;\dot S_i^h.

U_d = v; D_d = design direction        # tangent seeds
for scalar operation z = f(x1, ..., xn) in forward order:
    z_d = sum(partial(f, xk) * xk_d)
# z = x*y: z_d = y*x_d + x*y_d
# evaluate the same reconstruction, face, source and boundary operations
return R_d, J_d                          # residual and objective tangents

代码清单和播放器沿用 Tapenade 的命名：x_d 即 $\dot x$ ，x_b 即 $\overline x$ （生成的子程序称为 residual_d、residual_b）。播放器以固定的非零状态种子 $v$ 和参数种子 $\dot D=0$ 逐步展示 $\dot R_h=A_hv$ （表中的 R_d）的装配。第 10 节才求解设计引起的状态响应；一个算子调用与一次线性求解不是同一件事。

Code listings and players follow Tapenade's naming: x_d is $\dot x$ and x_b is $\overline x$ (the generated routines are residual_d and residual_b). The player uses a fixed nonzero state seed $v$ and the parameter seed $\dot D=0$ to display, step by step, the assembly of $\dot R_h=A_hv$ (R_d in the table). Section 10 solves the state response to a design direction; an operator application is not a linear solve.

算例：一个内部面的前向传播

Worked example: one interior face, forward

取页顶默认参数（N = 12，二阶，加热，常物性，温度跟踪）的收敛解，看第 $f=6$ 个面（ $x_f=0.4583$ ，左右是节点 5、6），设计方向 $\dot D=e_2$ （只改变 $\mu_{\rm ref}$ ），状态种子 $v$ 取第 10 节对这个方向解出的状态 tangent。下表每个量给出基态值，以及 tangent 的两部分：固定状态只改设计（ $v=0$ ，属于 $C_h\dot D$ ），与固定设计只改状态（ $\dot D=0$ ，属于 $A_hv$ ）。

Take the converged solution at the page defaults (N = 12, second order, heated, constant transport, tracking) and face $f=6$ ( $x_f=0.4583$ , between nodes 5 and 6), the design direction $\dot D=e_2$ (only $\mu_{\rm ref}$ changes), and as the state seed $v$ the state tangent that Section 10 solves for this direction. For each quantity the table gives the base value and the two parts of its tangent: the design varied with the state fixed ( $v=0$ , part of $C_h\dot D$ ) and the state varied with the design fixed ( $\dot D=0$ , part of $A_hv$ ).

量quantity	基态值base value	来自 $\dot D=e_2$ （ $v=0$ ）from $\dot D=e_2$ ( $v=0$ )	来自 $v$ （ $\dot D=0$ ）from $v$ ( $\dot D=0$ )	tangent 合计tangent, total
$G_f^u$	4.042×10⁻³	0	−0.01832	−0.01832
$G_f^T$	8.837×10⁻³	0	−0.04383	−0.04383
$u_f^v$	0.3512	0	−5.108×10⁻³	−5.108×10⁻³
$T_f^v$	1.005	0	−0.01791	−0.01791
$\mu_f$	0.1500	1.000	0	1.000
$\kappa_f$	0.7292	4.861	0	4.861
$\tau_f$	8.084×10⁻⁴	5.389×10⁻³	−3.664×10⁻³	1.725×10⁻³
$q_f$	−6.443×10⁻³	−0.04296	0.03196	−0.01099
$F^v_{f,3}$	6.727×10⁻³	0.04485	−0.03325	0.01159
$\widehat F^c_{f,3}$	1.252	0	−0.01795	−0.01795
$H_{f,2}$	1.123	−5.389×10⁻³	2.925×10⁻³	−2.464×10⁻³
$H_{f,3}$	1.245	−0.04485	0.01531	−0.02954

读法：常物性时 $\mu_f=\mu_{\rm ref}$ ，所以 $\dot\mu_f=\dot\mu_{\rm ref}=1$ ， $\dot\kappa_f=c_p/Pr=4.861$ ；固定状态时梯度不变， $\dot\tau_f=\frac43G_f^u=5.389\times10^{-3}$ ， $\dot q_f=-G_f^T\dot\kappa_f=-0.04296$ 。状态的响应使速度与温度梯度都变小（黏度与导热增大，剖面变平），抵消了一部分：合计 $\dot\tau_f=1.725\times10^{-3}$ 、 $\dot q_f=-0.01099$ 。对流通量不含设计，它的 tangent 全部来自 $v$ 。

Reading it: with constant transport $\mu_f=\mu_{\rm ref}$ , so $\dot\mu_f=\dot\mu_{\rm ref}=1$ and $\dot\kappa_f=c_p/Pr=4.861$ ; with the state fixed the gradients do not move, $\dot\tau_f=\frac43G_f^u=5.389\times10^{-3}$ and $\dot q_f=-G_f^T\dot\kappa_f=-0.04296$ . The response of the state makes both the velocity and the temperature gradient smaller (more viscosity and conduction flatten the profiles) and cancels part of this: in total $\dot\tau_f=1.725\times10^{-3}$ and $\dot q_f=-0.01099$ . The convective flux contains no design; its tangent comes entirely from $v$ .

把所有面的 $\dot H_f$ 分发到两侧、再减去源项的 tangent，就得到 $\dot R_h=A_hv+C_h\dot D$ ；因为 $v$ 是解出的 tangent，结果在线性求解的精度内为零。

Scattering the $\dot H_f$ of every face to its two sides and subtracting the source tangents gives $\dot R_h=A_hv+C_h\dot D$ ; since $v$ is the solved tangent, the result is zero to the accuracy of the linear solve.

图 7 用解出的 tangent 逐面播放前向扫描：设计方向 $\dot D=e_k$ 取页顶的“前向／连续对照方向”， $v$ 由第 10 节的求解器解出，按钮与本页其它播放器相同。每一步高亮当前面读取的存储状态：二阶对流重构读 $L-1$ 到 $R+1$ ，黏性部分只读 $L$ 、 $R$ （端面读三个节点）；然后把 $\dot H_f$ 加到左边控制体、从右边减去。条形是累加中的 $\dot R_h$ ：一个控制体收到两侧面的贡献之后就回到零附近，这正是 tangent 方程 $A_hv+C_h\dot D=0$ 逐行成立；热源方向要等最后减去源项的 tangent 才归零，格点页的五行最后被约束的 tangent 覆盖。表中给出当前面的 $\dot H_f$ 及其来自 $\dot D$ 与来自 $v$ 的两部分，最后一步给出 $\dot J_h=j_U^Tv+j_D^T\dot D=(g_h)_k$ 。

Figure 7 plays the forward sweep face by face with the solved tangent: the design direction $\dot D=e_k$ is the "tangent / reference direction" chosen at the top of the page, $v$ comes from the solver of Section 10, and the buttons are those of the page's other players. Each step highlights the stored states the current face reads: the second-order convective reconstruction reads $L-1$ to $R+1$ , the viscous part only $L$ and $R$ (an end face reads three nodes); then $\dot H_f$ is added to the left control volume and subtracted from the right one. The bars are $\dot R_h$ as it accumulates: once a control volume has received both of its faces it returns to near zero, which is the tangent equation $A_hv+C_h\dot D=0$ holding row by row; a heat-source direction returns to zero only after the source tangents are subtracted, and on the node page the five rows are finally overwritten by the constraint tangents. The table gives the current face's $\dot H_f$ and its parts from $\dot D$ and from $v$ ; the last step gives $\dot J_h=j_U^Tv+j_D^T\dot D=(g_h)_k$ .

图 7 · 前向扫描：用解出的 tangent

v=\dot{\mathbf U}_h

与

\dot D=e_k

逐面装配

\dot R_h

（在本页实时计算，跟随页顶的算例、目标与方向）。紫线是当前面，蓝框是对流重构读取的状态，绿色填充是黏性部分读取的状态，箭头是向两侧控制体的分发；条形按分量各自缩放。名称沿用 Tapenade：R_d 即

\dot R_h

，H_d 即

\dot H_f

，J_d 即

\dot J_h

。Figure 7 · The forward sweep:

\dot R_h

assembled face by face with the solved tangent

v=\dot{\mathbf U}_h

and

\dot D=e_k

(computed live on this page for the case, objective and direction selected at the top). The violet line is the current face, blue outlines mark the states the convective reconstruction reads, green fills the states the viscous part reads, and the arrows are the scatter to the two control volumes; each component's bars have their own scale. The names follow Tapenade: R_d is

\dot R_h

, H_d is

\dot H_f

, J_d is

\dot J_h

第 9 节Section 9

伴随：先反转边界覆盖，再反转装配

Adjoint: undo the boundary overwrite, then reverse the assembly

本节实验：逐步播放反向累加，数一数同一存储状态从几个面、几条路径收到贡献；共享输入必须累加，不能覆盖。

Experiment: Step through the reverse accumulation and count how many faces and paths contribute to one stored state. Shared inputs must accumulate, not overwrite.

反向传播计算什么

What reverse propagation computes

\Phi=\overline R_h^TR_h+\overline J_hJ_h,\qquad \overline{\mathbf U}_h=A_h^T\overline R_h+\overline J_hj_U,\qquad\overline D=C_h^T\overline R_h+\overline J_hj_D.

残差种子 $\overline R_h$ 是任意的 $3M$ 维向量， $\overline J_h$ 是目标种子；上横线表示标量 $\Phi$ 对某个计算节点的偏导，即反向模式导数，不是连续伴随场。取 $(\overline R_h,\overline J_h)=(\overline R_h,0)$ 得到矩阵转置乘积 $\overline{\mathbf U}_h=A_h^T\overline R_h$ ，取 $(0,1)$ 得到目标梯度 $j_U,j_D$ 。第 11 节的伴随解 $\psi_h$ 就是一个特定的残差种子：取 $\overline R_h=\psi_h$ 、 $\overline J_h=0$ ，反向扫描给出 $\overline{\mathbf U}_h=A_h^T\psi_h=-j_U$ ；取 $\overline J_h=1$ ，则 $\overline{\mathbf U}_h=0$ ，设计伴随 $\overline D=C_h^T\psi_h+j_D=g_h$ 就是总设计梯度。

The residual seed $\overline R_h$ is an arbitrary vector of length $3M$ and $\overline J_h$ seeds the objective. An overbar is a partial derivative of the scalar $\Phi$ with respect to a computational node, a reverse-mode derivative, not a continuous adjoint field. Seeds $(\overline R_h,\overline J_h)=(\overline R_h,0)$ give transpose products $\overline{\mathbf U}_h=A_h^T\overline R_h$ ; $(0,1)$ gives the objective gradients $j_U,j_D$ . The adjoint solution $\psi_h$ of Section 11 is one particular residual seed: with $\overline R_h=\psi_h$ and $\overline J_h=0$ the reverse sweep gives $\overline{\mathbf U}_h=A_h^T\psi_h=-j_U$ ; with $\overline J_h=1$ it gives $\overline{\mathbf U}_h=0$ , and the design adjoint $\overline D=C_h^T\psi_h+j_D=g_h$ is the total design gradient.

\overline R_h^TA_hv=\sum_r\sum_j\overline R_{h,r}\frac{\partial R_{h,r}}{\partial U_{h,j}}v_j=\sum_jv_j\sum_r\frac{\partial R_{h,r}}{\partial U_{h,j}}\overline R_{h,r}=v^TA_h^T\overline R_h.

转置来自有限求和换序，没有分部积分。参数方向同样换序，就得到前向与反向模式之间的对偶（点积）恒等式：

The transpose follows from exchanging finite sums, without integration by parts. Exchanging the sums for the parameter direction as well gives the duality (dot-product) identity between forward and reverse mode:

\overline R_h^T\dot R_h+\overline J_h\dot J_h=\overline{\mathbf U}_h^Tv+\overline D^T\dot D.

它对任意种子 $v,\dot D,\overline R_h,\overline J_h$ 都成立，是检查转置实现的点积测试。反向时对每个标量操作倒序应用链式法则：先把全部累加器清零（ $\overline x\leftarrow0$ ），再逐条执行 $\overline x\;{+}\!{=}\;\ldots$ ；共享输入必须累加。前向中一项被多次读取，反向就收到多条贡献。

It holds for any seeds $v,\dot D,\overline R_h,\overline J_h$ and is the dot-product test of the transpose implementation. In reverse mode apply the chain rule in reverse operation order: clear every accumulator ( $\overline x\leftarrow0$ ), then execute $\overline x\;{+}\!{=}\;\ldots$ statement by statement; shared inputs accumulate. Multiple forward reads become multiple reverse contributions.

点积测试能发现什么、不能发现什么，可以直接做实验（图 8）。第 8 节说过，冻结耗散速度会得到另一套线性算子。只在 tangent 里冻结，adjoint 仍对精确残差反向，两侧立刻不再相等；把 tangent 与 adjoint 一起冻结，两侧又相等到舍入误差，测试通过，可这对算子求的是另一个函数的导数，梯度与复数步长相差 $1\times10^{-6}$ 到 $6\times10^{-5}$ （相对）。点积测试检查的是“互为转置”；导数是否完整，要靠第 13、14 节重新求解的差分与复数步长。

What the dot-product test can and cannot find is a direct experiment (Figure 8). Section 8 noted that freezing the dissipation speed gives another linear operator. Freeze it in the tangent only, keep the adjoint exact, and the two sides disagree at once; freeze it in both, and the two sides agree to round-off again, so the test passes, yet this pair differentiates another function and its gradient differs from the complex step by $1\times10^{-6}$ to $6\times10^{-5}$ (relative). The dot-product test checks transposition; whether the derivative is complete takes the re-solved differences and the complex step of Sections 13 and 14.

图 8 · 点积测试证明的是“互为转置”，不是“导数正确”（页顶默认参数，N = 12）。(a) 8 组随机种子

(v,\dot D,\overline R_h,\overline J_h)

下，对偶恒等式两侧的相对差。圆点：本页的精确 tangent 与 adjoint，差在

10^{-15}

左右，是舍入误差。叉：tangent 冻结耗散速度（删去

-\tfrac12\dot\alpha_f(U_f^+-U_f^-)

）而 adjoint 仍精确，两侧不再相等，差

4\times10^{-7}

到

4\times10^{-5}

，测试发现了错误。方块：tangent 与 adjoint 一起冻结，两侧又相等到舍入误差，测试通过。(b) 用 LU 求解两种伴随得到的六个梯度分量，与第 14 节复数步长的相对差：精确伴随最多

6.7\times10^{-14}

，一起冻结的伴随约

1\times10^{-6}

到

6\times10^{-5}

。一起冻结的一对算子互为转置，却是另一个函数的导数；只有第 13、14 节的差分或复数步长能发现。低于

10^{-17}

的值画在

10^{-17}

。

Figure 8 · The dot-product test proves transposition, not correctness (page defaults, N = 12). (a) The relative gap between the two sides of the duality identity for 8 random seed sets

(v,\dot D,\overline R_h,\overline J_h)

. Dots: this page's exact tangent and adjoint; the gap is about

10^{-15}

, round-off. Crosses: the tangent freezes the dissipation speed (it drops

-\tfrac12\dot\alpha_f(U_f^+-U_f^-)

) while the adjoint stays exact; the sides no longer agree, the gap is

4\times10^{-7}

4\times10^{-5}

, and the test finds the error. Squares: tangent and adjoint freeze it together; the sides agree to round-off again and the test passes. (b) The six gradient components from both adjoints, each solved by LU, against the complex step of Section 14: the exact adjoint at most

6.7\times10^{-14}

, the jointly frozen one about

1\times10^{-6}

6\times10^{-5}

. The jointly frozen pair are transposes of each other and yet the derivative of another function; only the differences of Sections 13 and 14 can tell. Values below

10^{-17}

are drawn at

10^{-17}

先反转边界覆盖，再反转装配

Reverse boundary overwrites before assembly

\begin{gathered}\overline R^{\rm raw}\;{+}\!{=}\;M_{\rm bc}^T\overline R_h,\qquad\overline\zeta\;{+}\!{=}\;E_{\rm bc}^T\overline R_h,\\\overline H_f\;{+}\!{=}\;\overline R^{\rm raw}_{f-1}-\overline R^{\rm raw}_f\quad(0<f<M),\qquad\overline S_i^h\;{+}\!{=}\;-\overline R^{\rm raw}_i.\end{gathered}

\overline H_0\;{+}\!{=}\;-\overline R^{\rm raw}_0,\quad\overline H_M\;{+}\!{=}\;\overline R^{\rm raw}_{M-1},\quad A_h^T\overline R_h=A_{\rm raw}^T\overline R^{\rm raw}+\Bigl(\frac{\partial\zeta}{\partial\mathbf U_h}\Bigr)^T\overline\zeta.

$\overline R^{\rm raw}$ 是原始平衡行的伴随：掩码已把被覆盖行的种子清零。边界面只有一侧接控制体，所以种子只有一项。强边界约束的反传单独累加到状态与参数；不能把已经覆盖掉的守恒行又反传一次。格心版掩码为单位阵且没有约束插入，但边界通量本身仍需反传。

$\overline R^{\rm raw}$ is the adjoint of the raw balance rows; the mask has already zeroed the seeds of overwritten rows. A boundary face touches only one control volume, so its seed has a single contribution. Reverse strong constraints separately into states and parameters; do not differentiate overwritten balances again. Cells use an identity mask and no constraint insertion, but still reverse their boundary fluxes.

内部黏性面的逐项反传

Reverse an interior viscous face term by term

令 $\eta=\overline H_f$ ，其三个分量按质量、动量、能量排列。因为净通量减去黏性通量，黏性部分先收到 $-\eta$ 。下面各累加器先清零（ $\overline\tau_f,\overline u_f^v,\overline q_f,\ldots\leftarrow0$ ），再按显示顺序执行，即前向运算的倒序：

Let $\eta=\overline H_f$ , ordered as mass, momentum and energy. Net flux subtracts viscosity, so the viscous path starts with $-\eta$ . Clear these accumulators ( $\overline\tau_f,\overline u_f^v,\overline q_f,\ldots\leftarrow0$ ) and execute them in the displayed order, the reverse of the forward operations:

\overline\tau_f\;{+}\!{=}\;-\eta_2-u_f^v\eta_3,\quad\overline u_f^v\;{+}\!{=}\;-\tau_f\eta_3,\quad\overline q_f\;{+}\!{=}\;\eta_3,

\overline G_f^u\;{+}\!{=}\;\frac43\mu_f\overline\tau_f,\quad\overline\mu_f\;{+}\!{=}\;\frac43G_f^u\overline\tau_f,\quad\overline G_f^T\;{+}\!{=}\;-\kappa_f\overline q_f,\quad\overline\kappa_f\;{+}\!{=}\;-G_f^T\overline q_f,

\overline\mu_f\;{+}\!{=}\;\frac{c_p}{Pr}\overline\kappa_f,\quad\overline{Pr}\;{+}\!{=}\;-\frac{\kappa_f}{Pr}\overline\kappa_f,\quad\overline T_f^v\;{+}\!{=}\;\mu_{T,f}\overline\mu_f,\quad\overline{\mu}_{\rm ref}\;{+}\!{=}\;\frac{\mu_f}{\mu_{\rm ref}}\overline\mu_f.

\overline U_f^v\;{+}\!{=}\;u_U(U_f^v)^T\overline u_f^v+T_U(U_f^v)^T\overline T_f^v,

\overline U_L\;{+}\!{=}\;\frac12\overline U_f^v-\frac{u_U(U_L)^T\overline G_f^u+T_U(U_L)^T\overline G_f^T}{h},\quad\overline U_R\;{+}\!{=}\;\frac12\overline U_f^v+\frac{u_U(U_R)^T\overline G_f^u+T_U(U_R)^T\overline G_f^T}{h}.

这几行展示了面平均和原始量差分的两条反向路径。常物性模式只令 $\mu_{T,f}=0$ ，不能删除对 $\mu_{\rm ref}$ 的设计导数。格点页的加热出口是约束行 $p_{M-1}-\tau_M-D_5$ 与 $q_M$ ：它们的种子先从被替换的守恒行取出（第 4 节），再按约束的导数反传。约束里的 $\tau_M$ 、 $q_M$ 在格点页仍是本构量（出口的三点梯度，物性取节点 $M-1$ 的温度），所以 $\overline\tau_M$ 、 $\overline q_M$ 接着沿上面的黏性链反传（见下面的边界段），只是没有内部面的面平均与 Rusanov 耗散。

These equations expose both reverse paths: face averaging and primitive differences. Constant transport sets only $\mu_{T,f}=0$ , retaining the design derivative with respect to $\mu_{\rm ref}$ . On the node page the heated outlet is the constraint rows $p_{M-1}-\tau_M-D_5$ and $q_M$ : their seeds are first taken from the replaced balance rows (Section 4) and then reversed through the constraint derivatives. The $\tau_M$ and $q_M$ inside the constraints are still constitutive on the node page (the three-point outlet gradient, the transport at the temperature of node $M-1$ ), so $\overline\tau_M$ and $\overline q_M$ then continue along the viscous chain above (see the boundary paragraph below), only without the face average and the Rusanov dissipation of an interior face.

对流与二阶重构的转置

Transpose convection and second-order reconstruction

\overline\alpha_f\;{+}\!{=}\;-\frac12\eta^T(U_f^+-U_f^-),

\overline U_f^-\;{+}\!{=}\;\frac12F_U^c(U_f^-)^T\eta+\frac{\alpha_f}{2}\eta+\frac12(u_U(U_f^-)+c_U(U_f^-))^T\overline\alpha_f,

\overline U_f^+\;{+}\!{=}\;\frac12F_U^c(U_f^+)^T\eta-\frac{\alpha_f}{2}\eta+\frac12(u_U(U_f^+)+c_U(U_f^+))^T\overline\alpha_f,\quad c_U=\frac{\gamma}{2c}T_U.

\overline U_j\;{+}\!{=}\;W_{Lj}(x_f)\overline U_f^-+W_{Rj}(x_f)\overline U_f^+.

黏性通量的转置（分块形式）

The viscous flux transposed (block form)

把第 8 节第四步的分块 Jacobian 转置，就是上面逐项反传的紧凑写法（ $\overline F_f^v=-\overline H_f$ ；常物性时第一块只剩 $\tau_fu_U^T\overline F^v_{f,3}$ ）：

Transposing the block Jacobians of Step 4 in Section 8 gives the term-by-term reversal above in compact form ( $\overline F_f^v=-\overline H_f$ ; with constant transport the first block reduces to $\tau_fu_U^T\overline F^v_{f,3}$ ):

\overline U_f^v\;{+}\!{=}\;\Bigl(\frac{\partial F^v}{\partial U^v}\Bigr)^T\overline F_f^v,\quad\overline G_f^u\;{+}\!{=}\;\tfrac43\mu_f\bigl(\overline F^v_{f,2}+u_f^v\overline F^v_{f,3}\bigr),\quad\overline G_f^T\;{+}\!{=}\;\kappa_f\overline F^v_{f,3},\quad\overline\mu_{\rm ref}\;{+}\!{=}\;\frac{(F_f^v)^T\overline F_f^v}{\mu_{\rm ref}},\quad\overline{Pr}\;{+}\!{=}\;\frac{q_f}{Pr}\overline F^v_{f,3}.

梯度模板的转置

The transpose of the gradient stencil

第四步的模板 $G_f^\theta=\sum_jw_{fj}\theta(U_j)$ 反过来就是分发：每个模板点按同一个权重收到一份（图 9）：

Reversed, the stencil $G_f^\theta=\sum_jw_{fj}\theta(U_j)$ of Step 4 scatters: every stencil point receives a share with the same weight (Figure 9):

\overline U_j\;{+}\!{=}\;w_{fj}\,\theta_U(U_j)^T\,\overline G_f^\theta\quad(j\text{ in the stencil of }f,\ \theta=u,T).

内部面 $w=\mp1/h$ ，所以 $\overline U_L$ 减、 $\overline U_R$ 加，正是上面最后一行的两个 $1/h$ 项。一个存储状态出现在几个面的模板里，就从几个面收到贡献；这就是共享输入必须累加的原因。

At an interior face $w=\mp1/h$ , so $\overline U_L$ decreases and $\overline U_R$ increases: these are the two $1/h$ terms of the last line above. A stored state that appears in the stencils of several faces receives a contribution from each; this is why shared inputs must accumulate.

图 9 · 面梯度是存储状态的有限模板，tangent 沿模板收集，adjoint 沿同一模板的转置分发（格点页，N = 12，

h=1/12

，权重

w_j

按实际数值标出）。左：前向

\dot G_f^\theta=\sum_jw_j\,\theta_U(U_j)v_j

。右：反向

\overline U_j\;{+}\!{=}\;w_j\,\theta_U(U_j)^T\,\overline G_f^\theta

。同一组权重、方向相反：这就是“伴随是模板的转置”。

\theta

取

u

或

T

；完整的面依赖集合见图 B.1。

Figure 9 · A face gradient is a finite stencil of stored states: the tangent gathers along it and the adjoint scatters along its transpose (node page, N = 12,

h=1/12

, the weights

w_j

printed with their values). Left: forward,

\dot G_f^\theta=\sum_jw_j\,\theta_U(U_j)v_j

. Right: reverse,

\overline U_j\;{+}\!{=}\;w_j\,\theta_U(U_j)^T\,\overline G_f^\theta

. The same weights in the opposite direction: this is what "the adjoint is the transpose of the stencil" means.

\theta

u

T

; Figure B.1 gives the complete face dependencies.

边界面、边界行、源项与目标的反传

Boundary faces, boundary rows, sources and the objective in reverse

格点页先处理五个约束行：它们的种子就是 $\overline R_h$ 在这五行的分量，随后这五行的原始平衡被掩码清零（本节开头）。按第 8 节第五步倒序：

The node page first handles the five constraint rows: their seeds are the components of $\overline R_h$ on these rows, and the raw balances of the rows are then masked (the start of this section). In the reverse order of Step 5 in Section 8:

\overline\zeta_{\rm in}=(\overline R_{h,0},\overline R_{h,1},\overline R_{h,2})^T:\qquad\overline U_0\;{+}\!{=}\;\overline\zeta_{{\rm in},1}e_1+\overline\zeta_{{\rm in},2}\,u_U(U_0)^T+\overline\zeta_{{\rm in},3}\,T_U(U_0)^T,\qquad\overline D_4\;{+}\!{=}\;-\overline\zeta_{{\rm in},3},

\overline\zeta_{\rm out}=(\overline R_{h,3M-2},\overline R_{h,3M-1})^T:\qquad\overline U_{M-1}\;{+}\!{=}\;\overline\zeta_{{\rm out},1}\,p_U(U_{M-1})^T,\quad\overline\tau_M\;{+}\!{=}\;-\overline\zeta_{{\rm out},1},\quad\overline q_M\;{+}\!{=}\;\overline\zeta_{{\rm out},2},\quad\overline D_5\;{+}\!{=}\;-\overline\zeta_{{\rm out},1}\ \text{(heat)}.

$\overline\tau_M$ 、 $\overline q_M$ 再按上面的黏性链反传（物性在节点 $M-1$ 的温度上求值），梯度沿三点单侧模板分发到 $M-3,M-2,M-1$ 。温度跟踪目标 $j_D=0$ ，所以 $(g_h)_4=-\psi_{h,2}$ 、 $(g_h)_5=-\psi_{h,3M-2}$ ：入口温度约束行与出口应力约束行的乘子取负号，就是这两个边界数据的梯度。

$\overline\tau_M$ and $\overline q_M$ then reverse through the viscous chain above (with the transport evaluated at the temperature of node $M-1$ ), and the gradients scatter along the three-point one-sided stencil to $M-3,M-2,M-1$ . The tracking objective has $j_D=0$ , so $(g_h)_4=-\psi_{h,2}$ and $(g_h)_5=-\psi_{h,3M-2}$ : the multipliers of the inlet-temperature and outlet-stress constraint rows, negated, are the gradients with respect to these two boundary data.

源项与目标在前向程序里排在面循环之后，所以反向时先于各个面反传（仅在五个约束行之后）：

Sources and the objective come after the face loop in the forward program, so in reverse they are reversed before the faces (after the five constraint rows only):

\overline S_i^h\;{+}\!{=}\;-\overline R^{\rm raw}_i,\qquad\overline D_k\;{+}\!{=}\;Q_{ik}\,\overline S_{i,3}^h\ \ (k=0,1),\qquad\overline{\mathbf U}_h\;{+}\!{=}\;\overline J_h\,j_U,\qquad\overline D\;{+}\!{=}\;\overline J_h\,j_D.

$j_U$ 、 $j_D$ 是目标自己的反向扫描，第 11 节列出三种目标的公式。

$j_U$ and $j_D$ are the objective's own reverse sweep; Section 11 lists them for the three objectives.

整个面循环：语句、tangent 与 adjoint

The whole face loop: statement, tangent and adjoint

下表把一个内部面按代码的执行顺序逐条列出（core.js 的 evaluate）。tangent 一列从上往下执行；adjoint 一列从下往上执行，每条都是累加，所有系数都在基态求值，每行的 $u_U,T_U,c_U,F^c_U$ 在该行写明的状态上求值， $\theta$ 取 $u$ 与 $T$ 。每个 adjoint 行都是同一行 tangent 的转置：对任意 $\dot x$ 、 $\overline y$ ， $\overline y^T\dot y=\overline x^T\dot x$ 。按顺序执行完，得到的就是这个面对 $A_h^T\overline R_h$ 与 $C_h^T\overline R_h$ 的贡献。边界面用第 8 节第五步与上面的边界反传代替第 1–4 行，并且没有 Rusanov 耗散（第 9 行不存在，第 10 行只剩边界状态的物理通量 $F^c(U)$ ）。格点页的加热出口照常用出口处的三点梯度执行第 5–8 行；出口节点的动量行与能量行随后被约束 $p_{M-1}-\tau_M-D_5$ 与 $q_M$ 替换，所以 $\tau_M$ 、 $q_M$ 只通过这两条约束进入梯度。

The table lists one interior face statement by statement in the order the code executes it (evaluate() in core.js). Execute the tangent column top to bottom and the adjoint column bottom to top; every adjoint statement accumulates, every coefficient is evaluated at the base state, the $u_U,T_U,c_U,F^c_U$ of a row at the state named in that row, and $\theta$ stands for $u$ and $T$ . Each adjoint row is the transpose of the tangent in the same row: $\overline y^T\dot y=\overline x^T\dot x$ for any $\dot x$ , $\overline y$ . Executed in order, they give this face's contribution to $A_h^T\overline R_h$ and $C_h^T\overline R_h$ . Boundary faces replace rows 1–4 by Step 5 of Section 8 and the boundary reversal above, and have no Rusanov dissipation (row 9 is absent and row 10 reduces to the physical flux $F^c(U)$ of the boundary state). On the node page the heated outlet runs rows 5–8 as usual with the three-point outlet gradient; the outlet node's momentum and energy rows are then replaced by the constraints $p_{M-1}-\tau_M-D_5$ and $q_M$ , so $\tau_M$ and $q_M$ enter the gradient only through those two constraints.

#	语句（前向执行顺序）statement (forward order)	tangent：自上而下tangent: top to bottom	adjoint：自下而上adjoint: bottom to top
1	语句statement $U_f^-=\sum_jW_{Lj}(x_f)U_j$ $U_f^+=\sum_jW_{Rj}(x_f)U_j$	tangenttangent $\dot U_f^\pm=\sum_jW_{\cdot j}(x_f)v_j$	adjointadjoint $\overline U_j\;{+}\!{=}\;W_{Lj}(x_f)\overline U_f^-+W_{Rj}(x_f)\overline U_f^+$
2	语句statement $U_f^v=\tfrac12(U_L+U_R)$	tangenttangent $\dot U_f^v=\tfrac12(v_L+v_R)$	adjointadjoint $\overline U_L\;{+}\!{=}\;\tfrac12\overline U_f^v$ $\overline U_R\;{+}\!{=}\;\tfrac12\overline U_f^v$
3	语句statement $u_f^v=u(U_f^v)$ $T_f^v=T(U_f^v)$	tangenttangent $\dot u_f^v=u_U\dot U_f^v$ $\dot T_f^v=T_U\dot U_f^v$	adjointadjoint $\overline U_f^v\;{+}\!{=}\;u_U^T\overline u_f^v+T_U^T\overline T_f^v$
4	语句statement $G_f^\theta=\bigl(\theta(U_R)-\theta(U_L)\bigr)/h$	tangenttangent $\dot G_f^\theta=\bigl(\theta_U(U_R)v_R-\theta_U(U_L)v_L\bigr)/h$	adjointadjoint $\overline U_R\;{+}\!{=}\;\theta_U(U_R)^T\overline G_f^\theta/h$ $\overline U_L\;{-}\!{=}\;\theta_U(U_L)^T\overline G_f^\theta/h$
5	语句statement $\mu_f=\mu(T_f^v;\mu_{\rm ref})$	tangenttangent $\dot\mu_f=\mu_{T,f}\dot T_f^v+\frac{\mu_f}{\mu_{\rm ref}}\dot\mu_{\rm ref}$	adjointadjoint $\overline T_f^v\;{+}\!{=}\;\mu_{T,f}\overline\mu_f$ $\overline\mu_{\rm ref}\;{+}\!{=}\;\frac{\mu_f}{\mu_{\rm ref}}\overline\mu_f$
6	语句statement $\kappa_f=c_p\mu_f/Pr$	tangenttangent $\dot\kappa_f=\frac{c_p}{Pr}\dot\mu_f-\frac{\kappa_f}{Pr}\dot{Pr}$	adjointadjoint $\overline\mu_f\;{+}\!{=}\;\frac{c_p}{Pr}\overline\kappa_f$ $\overline{Pr}\;{+}\!{=}\;-\frac{\kappa_f}{Pr}\overline\kappa_f$
7	语句statement $\tau_f=\tfrac43\mu_fG_f^u$	tangenttangent $\dot\tau_f=\tfrac43\bigl(\mu_f\dot G_f^u+G_f^u\dot\mu_f\bigr)$	adjointadjoint $\overline G_f^u\;{+}\!{=}\;\tfrac43\mu_f\overline\tau_f$ $\overline\mu_f\;{+}\!{=}\;\tfrac43G_f^u\overline\tau_f$
8	语句statement $q_f=-\kappa_fG_f^T$	tangenttangent $\dot q_f=-\kappa_f\dot G_f^T-G_f^T\dot\kappa_f$	adjointadjoint $\overline G_f^T\;{+}\!{=}\;-\kappa_f\overline q_f$ $\overline\kappa_f\;{+}\!{=}\;-G_f^T\overline q_f$
9	语句statement $\alpha_f=\tfrac12\bigl(u_f^-+c_f^-+u_f^++c_f^+\bigr)$	tangenttangent $\dot\alpha_f=\tfrac12\sum_\pm(u_U+c_U)(U_f^\pm)\,\dot U_f^\pm$	adjointadjoint $\overline U_f^\pm\;{+}\!{=}\;\tfrac12(u_U+c_U)(U_f^\pm)^T\overline\alpha_f$
10	语句statement $\widehat F_f^c=\tfrac12\bigl(F^c(U_f^-)+F^c(U_f^+)\bigr)$ $\quad-\tfrac12\alpha_f(U_f^+-U_f^-)$	tangenttangent $\dot{\widehat F}_f^c=\tfrac12\sum_\pm F^c_U(U_f^\pm)\dot U_f^\pm$ $\quad-\tfrac12\dot\alpha_f(U_f^+-U_f^-)-\tfrac12\alpha_f(\dot U_f^+-\dot U_f^-)$	adjointadjoint $\overline U_f^-\;{+}\!{=}\;\tfrac12F^c_U(U_f^-)^T\overline{\widehat F}_f^c+\tfrac12\alpha_f\overline{\widehat F}_f^c$ $\overline U_f^+\;{+}\!{=}\;\tfrac12F^c_U(U_f^+)^T\overline{\widehat F}_f^c-\tfrac12\alpha_f\overline{\widehat F}_f^c$ $\overline\alpha_f\;{+}\!{=}\;-\tfrac12(U_f^+-U_f^-)^T\overline{\widehat F}_f^c$
11	语句statement $F_f^v=\bigl(0,\ \tau_f,\ u_f^v\tau_f-q_f\bigr)^T$	tangenttangent $\dot F_f^v=\bigl(0,\ \dot\tau_f,\ \tau_f\dot u_f^v+u_f^v\dot\tau_f-\dot q_f\bigr)^T$	adjointadjoint $\overline\tau_f\;{+}\!{=}\;\overline F^v_{f,2}+u_f^v\overline F^v_{f,3}$ $\overline u_f^v\;{+}\!{=}\;\tau_f\overline F^v_{f,3}$ $\overline q_f\;{+}\!{=}\;-\overline F^v_{f,3}$
12	语句statement $H_f=\widehat F_f^c-F_f^v$	tangenttangent $\dot H_f=\dot{\widehat F}_f^c-\dot F_f^v$	adjointadjoint $\overline{\widehat F}_f^c\;{+}\!{=}\;\overline H_f$ $\overline F_f^v\;{+}\!{=}\;-\overline H_f$
13	语句statement $R^{\rm raw}_{f-1}\;{+}\!{=}\;H_f$ $R^{\rm raw}_f\;{-}\!{=}\;H_f$	tangenttangent $\dot R^{\rm raw}_{f-1}\;{+}\!{=}\;\dot H_f$ $\dot R^{\rm raw}_f\;{-}\!{=}\;\dot H_f$	adjointadjoint $\overline H_f\;{+}\!{=}\;\overline R^{\rm raw}_{f-1}-\overline R^{\rm raw}_f$

x_b = 0 for every reverse accumulator
R_b = residual seed; J_b = objective seed
for operation z = f(inputs) in reverse order:
    for input x in inputs:
        x_b += partial(f, x) * z_b
# final boundary replacements are encountered first
# each shared input receives all contributions
return U_b, D_b

播放器以固定种子 $\overline R_h$ 、 $\overline J_h=0$ 逐步累加 $\overline{\mathbf U}_h=A_h^T\overline R_h$ （表中的 U_b）；它是任意种子的转置乘积，不是已经求出的伴随解。内部面的局部点积应与全局点积分别检查：局部正确不能排除全局分发或边界掩码错误。附录 B 可以选择面，独立查看实际局部 JVP/VJP 数值。

The player accumulates $\overline{\mathbf U}_h=A_h^T\overline R_h$ (U_b in the table) step by step for a fixed seed $\overline R_h$ with $\overline J_h=0$ : a transpose product for an arbitrary seed, not a solved adjoint. Check local and global dot products separately: local correctness does not exclude assembly or boundary-mask errors. Appendix B selects a face and displays actual local JVP/VJP values.

算例：同一个面的反向传播

Worked example: the same face, backward

仍取第 8 节算例的第 $f=6$ 个面，种子取第 11 节解出的伴随： $\overline R_h=\psi_h$ 、 $\overline J_h=1$ 。两侧节点 5、6 的行都是内部行，没有被掩码，所以 $\overline H_f=\overline R^{\rm raw}_{f-1}-\overline R^{\rm raw}_f$ 就是这两个节点的乘子块之差。按上表从下往上执行（常物性， $\mu_{T,f}=0$ ）：

Keep face $f=6$ of the example in Section 8 and seed with the adjoint solved in Section 11: $\overline R_h=\psi_h$ , $\overline J_h=1$ . The rows of nodes 5 and 6 on either side are interior and unmasked, so $\overline H_f=\overline R^{\rm raw}_{f-1}-\overline R^{\rm raw}_f$ is the difference of their multiplier blocks. Execute the table above from the bottom up (constant transport, $\mu_{T,f}=0$ ):

量quantity	值value
$\overline H_{f,1}$	5.206×10⁻⁴
$\overline H_{f,2}$	1.003×10⁻⁴
$\overline H_{f,3}$	−2.325×10⁻⁴
$\overline\tau_f=-\overline H_{f,2}-u_f^v\overline H_{f,3}$	−1.863×10⁻⁵
$\overline u_f^v=-\tau_f\overline H_{f,3}$	1.880×10⁻⁷
$\overline q_f=\overline H_{f,3}$	−2.325×10⁻⁴
$\overline G_f^u=\tfrac43\mu_f\overline\tau_f$	−3.727×10⁻⁶
$\overline G_f^T=-\kappa_f\overline q_f$	1.695×10⁻⁴
$\overline\kappa_f=-G_f^T\overline q_f$	2.055×10⁻⁶
$\overline\mu_f=\tfrac43G_f^u\overline\tau_f+\tfrac{c_p}{Pr}\overline\kappa_f$	9.888×10⁻⁶
$\overline\mu_{\rm ref}\;{+}\!{=}\;\tfrac{\mu_f}{\mu_{\rm ref}}\overline\mu_f$	9.888×10⁻⁶
$\overline{Pr}\;{+}\!{=}\;-\tfrac{\kappa_f}{Pr}\overline\kappa_f$	−2.081×10⁻⁶

这个面对总梯度的贡献是 $\overline\mu_{\rm ref}\;{+}\!{=}\;9.888\times10^{-6}$ 与 $\overline{Pr}\;{+}\!{=}\;-2.081\times10^{-6}$ ，正好等于第 12 节的逐面公式 $-(F_f^v)^T\overline H_f/\mu_{\rm ref}$ 与 $-q_f\overline H_{f,3}/Pr$ 。局部对偶恒等式把两次传播连起来：用第 8 节算例的 $\dot H_f$ （两部分之和），

This face contributes $\overline\mu_{\rm ref}\;{+}\!{=}\;9.888\times10^{-6}$ and $\overline{Pr}\;{+}\!{=}\;-2.081\times10^{-6}$ to the total gradient, exactly the face formulas $-(F_f^v)^T\overline H_f/\mu_{\rm ref}$ and $-q_f\overline H_{f,3}/Pr$ of Section 12. The local duality identity ties the two propagations together: with the $\dot H_f$ of the example in Section 8 (both parts),

\overline H_f^T\dot H_f=7.284434454\times10^{-6},\qquad\sum_j\overline U_j^Tv_j+\overline\mu_{\rm ref}\,\dot\mu_{\rm ref}=7.284434454\times10^{-6},

其中 $\overline U_j$ 是这个面分发给节点 4–7 的贡献（对流重构读四个状态）。两边的相对差是 9×10⁻¹⁶，在舍入误差以内。对所有面、源项与目标求和，就是本节开头的全局恒等式。

where $\overline U_j$ is what this face scatters to nodes 4–7 (the convective reconstruction reads four states). The relative difference of the two sides is 9×10⁻¹⁶, within round-off. Summed over all faces, sources and the objective this is the global identity at the start of this section.

图 10 把整个反向扫描逐步播放出来，种子同样取 $\overline R_h=\psi_h$ 、 $\overline J_h=1$ （在本页实时计算，跟随页顶的算例与目标）。顺序与程序相反：先是被替换的五行（它们的种子交给约束，原始平衡行随后被掩码清零），再是目标与源项，然后面从 $f=M$ 倒退到 $0$ 。每一步把 $\overline H_f$ 沿当前面的模板分发到存储状态，并把设计伴随累加到六个分量上。条形是 $\overline{\mathbf U}_h$ ：一个状态收到所有读取它的面之后就回到零附近，这正是伴随方程 $A_h^T\psi_h+j_U=0$ 逐行成立；下方六个条形是 $\overline D$ 占最终值 $g_h$ 的比例：两个热源幅度在源项那一步得到， $\mu_{\rm ref}$ 与 $Pr$ 一个面一个面地累积， $D_4$ 与 $D_5$ 在约束那一步一次得到（三种目标都不显式含有它们）。

Figure 10 plays the whole reverse sweep step by step, again with the seeds $\overline R_h=\psi_h$ , $\overline J_h=1$ (computed live on this page for the case and objective selected at the top). The order is the reverse of the program: first the five replaced rows (their seeds go to the constraints, and the raw balance rows are then masked), then the objective and the sources, then the faces from $f=M$ back to $0$ . Each step scatters $\overline H_f$ along the current face's stencil to the stored states and adds its design adjoint to the six components. The bars are $\overline{\mathbf U}_h$ : a state returns to near zero once every face that reads it has been reversed, which is the adjoint equation $A_h^T\psi_h+j_U=0$ holding row by row; the six bars below are $\overline D$ as a fraction of its final value $g_h$ : the two heat amplitudes arrive at the sources step, $\mu_{\rm ref}$ and $Pr$ build up face by face, and $D_4$ and $D_5$ arrive in one step at the constraints (none of the three objectives has an explicit term in them).

图 10 · 反向扫描：以

\overline R_h=\psi_h

、

\overline J_h=1

按程序的逆序累加

\overline{\mathbf U}_h

与

\overline D

（在本页实时计算，跟随页顶的算例与目标）。紫线是当前面，红色箭头是分发，蓝框是经对流重构收到贡献的存储状态，绿色填充是经黏性模板收到贡献的存储状态；上方条形按分量各自缩放，下方条形是

\overline D_k/(g_h)_k

（只画 −0.25 到 1.25，右边的数字是

\overline D_k

）。名称沿用 Tapenade：U_b 即

\overline{\mathbf U}_h

，D_b 即

\overline D

，H_b 即

\overline H_f

。Figure 10 · The reverse sweep:

\overline{\mathbf U}_h

and

\overline D

accumulated in the reverse order of the program with

\overline R_h=\psi_h

\overline J_h=1

(computed live on this page for the case and objective selected at the top). The violet line is the current face, red arrows are the scatter, blue outlines mark the stored states that receive a contribution through the convective reconstruction and green fills those that receive one through the viscous stencil; the upper bars have one scale per component, the lower bars are

\overline D_k/(g_h)_k

(drawn only between −0.25 and 1.25; the number at the right is

\overline D_k

). The names follow Tapenade: U_b is

\overline{\mathbf U}_h

, D_b is

\overline D

, H_b is

\overline H_f

第 10 节Section 10

前向求解：一个设计方向一次线性求解

Tangent solve: one linear solve per design direction

从收敛方程得到设计 tangent

Obtain the design tangent from the converged equations

假定当前光滑分支已收敛且 $A_h$ 非奇异。沿参数方向 $\dot D$ 改变设计，原始解也随之改变。这里要找的是 $v=\dot{\mathbf U}_h=d\mathbf U_h(D+\epsilon\dot D)/d\epsilon|_0$ ，不是任选的向量。对每个零残差分量求导：

Assume the current smooth branch is converged and $A_h$ is nonsingular. Changing design in direction $\dot D$ changes the primal solution. We seek $v=\dot{\mathbf U}_h=d\mathbf U_h(D+\epsilon\dot D)/d\epsilon|_0$ , not an arbitrarily chosen vector. Differentiate each zero-residual component:

R_h(\mathbf U_h(D),D)=0\quad\Longrightarrow\quad\sum_j\frac{\partial R_{h,r}}{\partial U_{h,j}}v_j+\sum_{k=0}^{5}\frac{\partial R_{h,r}}{\partial D_k}\dot D_k=0,\quad A_hv=-C_h\dot D.

求一个参数的导数取 $\dot D=e_k$ ；同时改变多个参数则取一般方向。先以零状态种子求出右端，再解线性系统，最后把求出的状态方向和同一个参数方向一起传入目标求导。目标的显式参数项不能遗漏。

Use $\dot D=e_k$ for one parameter derivative, or a general direction for simultaneous parameter changes. Form the right-hand side with a zero state seed, solve the system, and finally differentiate the objective using both the solved state direction and the same parameter direction. Retain explicit objective parameter terms.

\widehat J_h(D)=J_h(\mathbf U_h(D),D),\qquad\dot J_h=\left.\frac{d}{d\epsilon}\widehat J_h(D+\epsilon\dot D)\right|_{\epsilon=0}=j_U^Tv+j_D^T\dot D.

solve R_h(U,D)=0 to the primal tolerance
L = exact_linearization(U,D)
D_d = e_k                                   # or any design direction
rhs = -L.jvp(U_d=0, D_d).R_d                # -C_h D_d
v = solve_linear(A_times=L.jvp(U_d=., D_d=0).R_d, rhs=rhs, preconditioner=B_h)
check norm(A_times(v)-rhs)
J_d = L.jvp(U_d=v, D_d).J_d                 # j_U^T v + j_D^T D_d

该过程得到所有状态对一个设计方向的响应。六个独立方向原则上需要六次 tangent 求解，但可以复用同一基态线性化和预处理。下图按密度、动量密度、能量密度绘制解向量；它不是速度和温度响应，若需要原始量响应，必须再用第 8 节的转换 Jacobian。

This computes the response of every state entry to one design direction. Six independent directions generally require six tangent solves, while sharing the base linearization and preconditioner. The plot shows density, momentum-density and energy-density responses, not velocity and temperature responses; obtain those using Section 8's conversion Jacobian.

如果原始残差尚未收敛，求出的敏感度对应当前线性化方程，不一定等于精确稳态分支的导数。如果 $A_h$ 接近奇异，即使线性残差很小，解误差仍可能被放大。此时需同时检查原始收敛、线性收敛和差分验证，而不能只看曲线是否平滑。

Before primal convergence, the sensitivity belongs to the current linearized equation and need not equal the exact steady-branch derivative. Near singularity, a small linear residual can still correspond to an amplified solution error. Check primal convergence, linear convergence and finite differences together rather than judging only smooth-looking curves.

右端 $-C_h\dot D$ ：每个设计方向从哪里进入

The right-hand side $-C_h\dot D$ : where each design direction enters

由第 8 节第四、五步， $C_h$ 的六列可以直接读出（固定状态，只让一个设计分量变化）：

Steps 4 and 5 of Section 8 give the six columns of $C_h$ directly (the state fixed, one design component varied):

\begin{gathered}\frac{\partial R_i^{\rm raw}}{\partial D_k}=-e_3Q_{ik}\ (k=0,1),\qquad\frac{\partial R_i^{\rm raw}}{\partial\mu_{\rm ref}}=-\frac{F^v_{i+1}-F^v_i}{\mu_{\rm ref}},\qquad\frac{\partial R_i^{\rm raw}}{\partial Pr}=-\frac{q_{i+1}-q_i}{Pr}e_3,\\\frac{\partial\zeta_{\rm out}}{\partial\mu_{\rm ref}}=\Bigl(-\frac{\tau_M}{\mu_{\rm ref}},\ \frac{q_M}{\mu_{\rm ref}}\Bigr)^T\ \text{(heat)},\qquad\frac{\partial\zeta_{{\rm out},2}}{\partial Pr}=-\frac{q_M}{Pr}\ \ \text{(heat)},\\\frac{\partial\zeta_{{\rm in},3}}{\partial D_4}=-1,\qquad\frac{\partial\zeta_{{\rm out},1}}{\partial D_5}=-1\ \text{(heat)},\end{gathered}

其中被约束替换的五行取约束的偏导（其余为零）；端面的黏性通量所在的行都被替换或只剩质量行，所以求和实际只到内部面。于是热源方向只动能量行， $\mu_{\rm ref}$ 动每个内部控制体的动量与能量行（以及出口两个约束行）， $Pr$ 只动能量行（包括出口的热流约束行）， $D_4$ 与 $D_5$ 各只动一个约束行。出口热流约束 $q_M=-\kappa_MG_M^T$ 的导热系数含 $\mu_{\rm ref}$ 与 $Pr$ ，所以它对二者的偏导 $q_M/\mu_{\rm ref}$ 与 $-q_M/Pr$ 是同一类；加热算例的收敛解上 $q_M=0$ （到牛顿求解的精度），这两项随之为零，这一列里非零的是出口应力行的 $-\tau_M/\mu_{\rm ref}$ 。

where the five rows replaced by constraints take the partial derivatives of the constraints (the others vanish); the rows of the end faces' viscous fluxes are replaced or keep only their mass row, so the sums effectively run over interior faces. A heat direction therefore moves only energy rows, $\mu_{\rm ref}$ the momentum and energy rows of every interior control volume (and the two outlet constraint rows), $Pr$ only energy rows (the outlet heat-flux constraint row included), and $D_4$ and $D_5$ one constraint row each. The conductivity in the outlet heat-flux constraint $q_M=-\kappa_MG_M^T$ contains $\mu_{\rm ref}$ and $Pr$ , so its derivatives $q_M/\mu_{\rm ref}$ and $-q_M/Pr$ are of one kind; in the heated case $q_M=0$ at the converged solution (to the accuracy of the Newton solve), both vanish with it, and the nonzero entry of these columns is the $-\tau_M/\mu_{\rm ref}$ of the outlet-stress row.

可是解出 $A_hv=-C_h\dot D$ 之后，每个方向都改变了所有状态。图 11 用第 8 节的转换 Jacobian 把六个解出的 tangent 换成温度 $\dot T_i=T_U(U_i)v_i$ ，并给出每个方向的 $\dot J_h=(g_h)_k$ ；六次求解共 206 次 GMRES 迭代（每个方向 34、34、34、34、35、35 次）。

After solving $A_hv=-C_h\dot D$ , however, every direction changes every state. Figure 11 converts the six solved tangents to temperature, $\dot T_i=T_U(U_i)v_i$ , with the conversion Jacobian of Section 8 and gives $\dot J_h=(g_h)_k$ for each direction; the six solves take 206 GMRES iterations (34, 34, 34, 34, 35, 35 per direction).

图 11 · 六个设计方向

\dot D=e_k

的解出的状态 tangent，换成温度：

\dot T_i=T_U(U_i)v_i

（第 8 节的转换 Jacobian；页顶默认参数，节点值，N = 12）。最上面是温度跟踪目标的失配

T-T^*

：它在 11 个节点上为负，只在最后一个节点略为正，所以使温度升高的方向（两个热源幅度、

Pr

、入口温度、出口应力）使目标减小，梯度为负；

\mu_{\rm ref}

使温度降低，梯度为正。每格给出

\dot J_h=(g_h)_k

，与第 11 节一次伴随求解得到的六个分量相同。

Figure 11 · The solved state tangents of the six design directions

\dot D=e_k

, converted to temperature:

\dot T_i=T_U(U_i)v_i

(the conversion Jacobian of Section 8; page defaults, nodal values, N = 12). At the top, the mismatch

T-T^*

of the tracking objective: it is negative at 11 of the 13 nodes and slightly positive only at the last node, so directions that raise the temperature (the two heat amplitudes,

Pr

, the inlet temperature, the outlet stress) lower the objective and have negative gradients, while

\mu_{\rm ref}

lowers the temperature and has a positive one. Each panel gives

\dot J_h=(g_h)_k

, the same six components as the one adjoint solve of Section 11.

第 11 节Section 11

伴随求解：一次求解，全部六个参数

Adjoint solve: one solve for all six parameters

有限维 Lagrangian：不是连续积分形式

A finite-dimensional Lagrangian, not a continuous integral

\mathcal L_h(\mathbf U_h,D,\psi_h)=J_h(\mathbf U_h,D)+\psi_h^TR_h(\mathbf U_h,D),\qquad\psi_h\in\mathbb R^{3M}.

$\psi_h$ 的每个分量对应一条最终离散残差行。它是有限维乘子向量，不是定义在连续空间上的函数。格点边界行若已改为约束，相应乘子就对应那条约束。这里采用加号约定，后面的负号均由此推导。

Each component of $\psi_h$ multiplies one final discrete residual row. It is a finite-dimensional multiplier vector, not a continuous spatial function. A node boundary row replaced by a constraint has a multiplier for that constraint. The plus convention here determines the signs below.

把每一项展开，再消去状态响应

Expand each term before eliminating the state response

\dot{\mathcal L}_h=j_U^Tv+j_D^T\dot D+\psi_h^T(A_hv+C_h\dot D)=(j_U+A_h^T\psi_h)^Tv+(j_D+C_h^T\psi_h)^T\dot D.

在收敛解上 $R_h=0$ ，所以 Lagrangian 等于目标。即使乘子随设计改变，额外项 $\dot\psi_h^TR_h$ 仍为零。选择乘子使未知状态响应的系数消失，就不再需要为每个参数单独求 $v$ 。

At convergence $R_h=0$ , so the Lagrangian equals the objective. Even if multipliers vary with design, the extra term $\dot\psi_h^TR_h$ vanishes. Choose multipliers to cancel the coefficient of the unknown state response and avoid solving a separate $v$ for every parameter.

\boxed{A_h^T\psi_h=-j_U},\qquad\boxed{g_h=j_D+C_h^T\psi_h},\qquad (g_h)_k=\frac{\partial J_h}{\partial D_k}+\sum_r\psi_{h,r}\frac{\partial R_{h,r}}{\partial D_k}.

$g_h$ 是六维总设计梯度， $j_U$ 是 $3M$ 维目标状态梯度，二者不共用字母 g。一个标量目标只需一次伴随求解；换目标一般要换右端并重新求解。换参数方向不改变伴随解，但若改变实际设计并重算原始流场，则必须重新线性化。

The total design gradient $g_h$ has six components; the state objective gradient $j_U$ has $3M$ components. They no longer share the letter g. One scalar objective needs one adjoint solve; a different objective generally changes the right-hand side and requires another solve. Selecting another parameter direction leaves the adjoint unchanged, but changing the actual design and re-solving the primal requires relinearization.

j_U^Tv=-\psi_h^TA_hv=\psi_h^TC_h\dot D,\qquad\dot J_h=j_U^Tv+j_D^T\dot D=g_h^T\dot D.

这是 tangent 与 adjoint 对同一离散问题给出相同方向导数的代数证明。没有连续极限假设，也没有分部积分。如果改用减号定义离散 Lagrangian，乘子会整体换号；总梯度仍应一致。这个符号约定不等于连续场与离散乘子的逐点对应关系。

This algebra proves equality of tangent and adjoint directional derivatives for the same discrete problem. It uses neither a continuum limit nor integration by parts. A minus-sign discrete Lagrangian changes multiplier signs while preserving total gradients. That sign convention does not establish a pointwise correspondence between continuous fields and discrete multipliers.

jU = L.vjp(R_b=0, J_b=1).U_b
AT_times(R_b) = L.vjp(R_b, J_b=0).U_b
psi = solve_linear(AT_times, rhs=-jU, preconditioner=transpose_GS)
check norm(AT_times(psi)+jU)
g = L.vjp(R_b=psi, J_b=1).D_b              # the design adjoint = g_h

三种目标的 $j_U$ 与 $j_D$

$j_U$ and $j_D$ of the three objectives

伴随右端是 $-j_U$ ，总梯度还要加 $j_D$ 。温度跟踪的 $j_U$ 已在第 7 节写出， $j_D=0$ 。两个端点目标只读一个面上的量，所以 $j_U$ 只有几个块非零，而且可能直接依赖设计：

The adjoint right-hand side is $-j_U$ , and the total gradient adds $j_D$ . Section 7 wrote $j_U$ of temperature tracking, whose $j_D=0$ . The two endpoint objectives read quantities of one face only, so only a few blocks of $j_U$ are nonzero, and they may depend on the design directly:

\begin{gathered}J_{\rm out}=T(U_{M-1}):\quad(j_U)_{M-1}=T_U(U_{M-1})^T,\quad j_D=0;\\J_q=q_0=-\kappa_0G_0^T,\ G_0^T=\frac{-3T_0+4T_1-T_2}{2h}:\quad(j_U)_0=\Bigl(\frac{3\kappa_0}{2h}-G_0^T\kappa_{T,0}\Bigr)T_U(U_0)^T,\ (j_U)_1=-\frac{2\kappa_0}{h}T_U(U_1)^T,\ (j_U)_2=\frac{\kappa_0}{2h}T_U(U_2)^T,\\ j_D=\Bigl(0,\ 0,\ \frac{q_0}{\mu_{\rm ref}},\ -\frac{q_0}{Pr},\ 0,\ 0\Bigr)^T.\end{gathered}

格点页入口热流的两个显式项来自 $\kappa_0\propto\mu_{\rm ref}/Pr$ ；入口温度只经约束行进入，没有显式项。页顶默认参数下 $q_0=-7.232\times10^{-3}$ ， $j_D=(0,0,-0.04822,0.01004,0,0)$ 。制造解算例的出口温度目标也是 $T(U_{M-1})$ 。

On the node page the two explicit terms of the inlet heat flux come from $\kappa_0\propto\mu_{\rm ref}/Pr$ ; the inlet temperature enters only through its constraint row and has no explicit term. At the page defaults $q_0=-7.232\times10^{-3}$ and $j_D=(0,0,-0.04822,0.01004,0,0)$ . The outlet-temperature objective of the manufactured case is $T(U_{M-1})$ as well.

用真残差判断伴随求解

Judge the adjoint solve using the true residual

r_\psi=A_h^T\widetilde\psi_h+j_U,\qquad \widetilde g_h^T\dot D-\dot J_h=-r_\psi^Tv,\qquad |r_\psi^Tv|\le\|r_\psi\|_2\|v\|_2.

此处 $\widetilde\psi_h$ 是尚未完全收敛的迭代解， $\widetilde g_h=j_D+C_h^T\widetilde\psi_h$ 是用它算出的梯度， $\dot J_h$ 是精确 tangent 给出的方向导数；等式假定使用同一精确基态且 tangent 已精确求解。它说明伴随残差与敏感度幅度共同决定方向导数误差。下图的乘子分量是离散行权重，不能直接读成连续密度、速度或温度伴随场。

Here $\widetilde\psi_h$ is an incompletely converged iterate, $\widetilde g_h=j_D+C_h^T\widetilde\psi_h$ the gradient computed from it and $\dot J_h$ the directional derivative of the exact tangent. The identity assumes the same exact base state and an exactly solved tangent. Adjoint residual and sensitivity magnitude jointly control directional-derivative error. The plotted multiplier components are weights on discrete rows, not direct continuous density, velocity or temperature adjoint fields.

乘子的含义：在一个控制体里加一点源项

What the multipliers mean: a small source in one control volume

在第 $r$ 行加一个源项 $s$ ，即把 $R_{h,r}$ 换成 $R_{h,r}-s$ ，再重新求解。对 $s$ 求导，状态的导数 $v=d\mathbf U_h/ds$ 满足 $A_hv-e_r=0$ ，于是

Add a source $s$ to row $r$ , that is, replace $R_{h,r}$ by $R_{h,r}-s$ , and re-solve. Differentiated with respect to $s$ , the state derivative $v=d\mathbf U_h/ds$ satisfies $A_hv-e_r=0$ , so

\frac{dJ_h}{ds}=j_U^Tv=j_U^TA_h^{-1}e_r=-\psi_h^Te_r=-\psi_{h,r}.

所以 $-\psi_{i,k}$ 是“每单位时间在控制体 $i$ 加入一单位第 $k$ 个守恒量（质量、动量或能量），目标改变多少”：伴随场是目标对源项的敏感度。热源梯度 $(g_h)_k=-\sum_i\psi_{i,3}Q_{ik}$ 就是能量行的敏感度按加热形状加权求和。图 12 画出三种目标的伴随场，并在三个节点重新求解加了源项的非线性问题来验证这个含义。温度跟踪的 $\psi_{i,3}>0$ ：任何位置加热都使目标减小，中部最有效；出口温度的 $\psi_{i,3}<0$ ，越靠近出口加热，出口温度升得越多；入口热流的 $\psi_{i,3}$ 在入口附近最大，那里的热量更多地逆着流动导出入口。第 16 节说明动量行与能量行的乘子在加密下逼近连续乘子 $-\lambda_\Pi$ 、 $-\lambda_{\mathcal H}$ 。被约束替换的五行，乘子的含义是约束目标值的敏感度：把 $T_0=D_4$ 改成 $T_0=D_4+s$ ， $dJ_h/ds=-\psi_{h,2}$ ，正是 $(g_h)_4$ 。

So $-\psi_{i,k}$ is how much the objective changes when one unit of the $k$ -th conserved quantity (mass, momentum or energy) per unit time is added to control volume $i$ : the adjoint field is the sensitivity of the objective to sources. The heat-source gradient $(g_h)_k=-\sum_i\psi_{i,3}Q_{ik}$ is the energy-row sensitivity weighted by the heating shape. Figure 12 draws the adjoint fields of the three objectives and checks this meaning by re-solving the nonlinear problem with a source added in three nodes. For temperature tracking $\psi_{i,3}>0$ : heating anywhere lowers the objective, most effectively in the middle; for the outlet temperature $\psi_{i,3}<0$ , and the closer to the outlet the heat is added, the more the outlet temperature rises; for the inlet heat flux $\psi_{i,3}$ is largest near the inlet, where more of the heat is conducted out through the inlet against the flow. Section 16 shows that the momentum and energy multipliers approach the continuous multipliers $-\lambda_\Pi$ , $-\lambda_{\mathcal H}$ under refinement. The five rows replaced by constraints carry the sensitivity to the constraint's target value: changing $T_0=D_4$ to $T_0=D_4+s$ gives $dJ_h/ds=-\psi_{h,2}$ , which is $(g_h)_4$ .

图 12 · 三种目标的离散伴随场（页顶默认参数，节点，N = 12）：每行一个目标，每列一个分量。曲线是

\psi_{i,1},\psi_{i,2},\psi_{i,3}

（质量、动量、能量行）；方块是独立的验证：在第 2、6、10 个节点的该行加入源项

s

，重新求解非线性原始问题，用中心差分得到

-dJ_h/ds

，与

\psi_{i,k}

最多相差 2×10⁻⁷（相对）。所以

-\psi_{i,k}

就是“在控制体

i

每单位时间加入一单位第

k

个守恒量，目标改变多少”。格点页被约束替换的五行（节点 0 的三行、节点 M−1 的动量与能量行）是约束乘子，含义不同，没有画出。

Figure 12 · The discrete adjoint fields of the three objectives (page defaults, nodes, N = 12): one objective per row, one component per column. The curves are

\psi_{i,1},\psi_{i,2},\psi_{i,3}

(the mass, momentum and energy rows); the squares are an independent check: a source

s

is added to that row of nodes 2, 6 and 10, the nonlinear primal is re-solved and central differences give

-dJ_h/ds

, which agrees with

\psi_{i,k}

to 2×10⁻⁷ (relative) at most. So

-\psi_{i,k}

is how much the objective changes when one unit of the

k

-th conserved quantity per unit time is added to control volume

i

. The five rows the node page replaces by constraints (the three of node 0, the momentum and energy rows of node M−1) hold constraint multipliers with another meaning and are not drawn.

代价：tangent 按设计分量计，adjoint 按目标计

Cost: the tangent counts design components, the adjoint counts objectives

tangent 每个设计方向解一次 $A_hv=-C_h\dot D$ ，得到所有目标沿这个方向的导数；adjoint 每个目标解一次 $A_h^T\psi_h=-j_U$ ，得到这个目标对所有设计分量的导数。两种求解的每次迭代代价大致相同：两次矩阵自由乘积（一次做 Krylov 步，一次检查真残差；每次都沿同一张记录了 4247 个中间量的计算图前向或反向扫描一次，两者都把每个中间量访问一次）加一次 GS 扫描，收尾时再多一次乘积（转置的 GS 扫描每次还要重新分解各个 $3\times3$ 对角块，所以略贵）。例如温度跟踪的伴随求解，35 次迭代是 71 次乘积和 35 次 GS 扫描。图 13 用本页默认参数的实际迭代次数比较。

The tangent solves $A_hv=-C_h\dot D$ once per design direction and obtains the derivatives of every objective along it; the adjoint solves $A_h^T\psi_h=-j_U$ once per objective and obtains the derivatives with respect to every design component. An iteration costs about the same in both: two matrix-free products (one for the Krylov step, one for the true-residual check, each a forward or reverse sweep over the 4247 values recorded in the same computational graph, visiting every value once) plus one GS sweep, and one more product at the end (the transposed GS sweep also refactors each $3\times3$ diagonal block on every call, so it is slightly dearer). The 35 iterations of the adjoint solve for temperature tracking, for example, are 71 products and 35 GS sweeps. Figure 13 compares the actual iteration counts at the page defaults.

图 13 · tangent 与 adjoint 的代价（页顶默认参数，N = 12；单位是 GMRES 迭代次数，每次迭代做两次完整的矩阵自由乘积（一次 Krylov 步，一次真残差检查，各沿记录了 4247 个中间量的计算图前向或反向扫描一次）和一次块 GS 扫描，收尾时再多一次乘积；每次迭代的结构在两条路线里相同，所以比较迭代次数就是比较代价，转置的 GS 扫描略贵）。(a) 一个目标：要

k

个设计分量的导数，tangent 要

k

次线性求解，累计 34、68、102、136、171、206 次迭代；adjoint 一次求解 35 次迭代就给出全部六个分量。(b) 三个目标、全部六个分量：tangent 的六次求解与目标无关，仍是 206 次（再加一个目标只需它自己的

j_U

、

j_D

：一次反向扫描和六个内积）；adjoint 每个目标一次求解，共 105 次。设计分量多用 adjoint，目标多、设计分量少用 tangent。

Figure 13 · The cost of tangent and adjoint (page defaults, N = 12; the unit is a GMRES iteration, each one making two complete matrix-free products (one for the Krylov step, one for the true-residual check, each a forward or reverse sweep over the 4247 values recorded in the computational graph) and one block-GS sweep, plus one more product at the end; an iteration has the same structure on both routes, so comparing iteration counts compares costs, the transposed GS sweep being slightly dearer). (a) One objective:

k

design components need

k

tangent solves, 34, 68, 102, 136, 171, 206 iterations cumulatively; one adjoint solve of 35 iterations gives all six. (b) Three objectives, all six components: the six tangent solves do not depend on the objective, still 206 iterations (a further objective only needs its own

j_U

and

j_D

: one reverse sweep and six dot products); the adjoint needs one solve per objective, 105 iterations in all. Many design components favour the adjoint; many objectives and few design components favour the tangent.

第 12 节Section 12

设计导数：热源、物性与边界参数

Design derivatives: heating, transport and boundary parameters

六个参数经过哪些计算节点

Where the six parameters enter the graph

两个热源幅度进入能量源项； $\mu_{\rm ref}$ 同时改变应力、导热及黏性功； $Pr$ 改变导热；入口温度与出口参数通过第 4 节的闭合进入。参数偏导固定全部状态槽位，状态随后如何变化由 tangent 或 adjoint 处理，不能在局部偏导中重复加入。

Two heat amplitudes enter the energy source. The parameter $\mu_{\rm ref}$ changes stress, conduction and viscous work, while $Pr$ changes conduction. Inlet temperature and the outlet parameter enter Section 4's closures. Parameter partials hold all state slots fixed; tangent or adjoint accounts for the resulting state response, which must not be counted again inside local partials.

热源梯度可以逐行写出

Heat-source gradients can be written row by row

Q_{ik}=\sum_{\pm}\omega_{i\pm}\sin^2((k+1)\pi x_{i\pm}),\quad\frac{\partial R_i^{\rm raw}}{\partial D_k}=-e_3Q_{ik},\quad (g_h)_k=\overline D_k=-\sum_i(M_{\rm bc}^T\psi_h)_{i,3}Q_{ik},\quad k=0,1.

负号来自残差减去源项，不是人为指定的伴随符号。本页目标与约束对两个热源幅度均无显式依赖，因此公式没有额外项。格点入口和出口的能量行已被替换，对应源项种子先被掩码置零；格心则保留所有控制体能量平衡。

The minus sign comes from subtracting the source in the residual, not an imposed adjoint sign. Objectives and constraints have no explicit dependence on the two heat amplitudes, so no extra terms appear here. Replaced inlet and outlet energy rows on nodes have their source seeds masked to zero; cells retain every volume energy balance.

物性与边界参数使用完整反传

Use full reverse accumulation for transport and boundary parameters

g_h=\overline D=j_D+C_{\rm raw}^TM_{\rm bc}^T\psi_h+\zeta_D^TE_{\rm bc}^T\psi_h.

这就是第 9 节以种子 $\overline R_h=\psi_h$ 、 $\overline J_h=1$ 反向扫描时在六个参数上累加出的设计伴随 $\overline D$ 。格心版最后一项为空，但边界通量对参数的影响已在中间项中。格点版入口温度约束对 $D_4$ 的显式偏导为 −1；出口应力或温度约束对 $D_5$ 的显式偏导同样为 −1。除此之外，还要保留物性和目标的直接参数路径。

This is the design adjoint $\overline D$ that the reverse sweep of Section 9 accumulates on the six parameters for the seeds $\overline R_h=\psi_h$ and $\overline J_h=1$ . The last term is empty for cells, but their boundary-flux parameter dependence is already in the middle term. For nodes, the inlet temperature constraint has explicit derivative −1 with respect to $D_4$ , and the outlet stress or temperature constraint has derivative −1 with respect to $D_5$ . Retain transport and direct objective paths as well.

逐面写出的设计导数

The design derivatives written face by face

把第 8 节第四、五步的偏导代入 $g_h=j_D+C_h^T\psi_h$ ，每个分量都能写成对面、控制体或约束行的求和，其中 $\overline H_f=\overline R^{\rm raw}_{f-1}-\overline R^{\rm raw}_f$ 由 $\psi_h$ 经第 9 节的装配转置得到：

Substituting the partial derivatives of Steps 4 and 5 of Section 8 into $g_h=j_D+C_h^T\psi_h$ writes each component as a sum over faces, control volumes or constraint rows, with $\overline H_f=\overline R^{\rm raw}_{f-1}-\overline R^{\rm raw}_f$ obtained from $\psi_h$ by the transposed assembly of Section 9:

\begin{gathered}(g_h)_k=-\sum_iQ_{ik}\overline R^{\rm raw}_{i,3}\ (k=0,1),\qquad(g_h)_2=j_{D,2}-\frac1{\mu_{\rm ref}}\Bigl(\sum_{0<f<M}\overline H_f^TF_f^v+\psi_{h,3M-2}\tau_M-\psi_{h,3M-1}q_M\Bigr),\\(g_h)_3=j_{D,3}-\frac1{Pr}\Bigl(\sum_{0<f<M}\overline H_{f,3}\,q_f+\psi_{h,3M-1}q_M\Bigr),\qquad(g_h)_4=j_{D,4}-\psi_{h,2},\qquad(g_h)_5=j_{D,5}-\psi_{h,3M-2}.\end{gathered}

这是加热算例的写法（制造解算例的出口约束 $u-0.35$ 、 $T-D_5$ 不含 $\mu_{\rm ref}$ 与 $Pr$ ，含 $\tau_M,q_M$ 的两项消失， $D_5$ 改由温度约束行的 $\psi_{h,3M-1}$ 给出）。 $\overline R^{\rm raw}=M_{\rm bc}^T\psi_h$ 已把五个约束行清零，两个端面因此不出现； $q_M$ 在收敛解上为零（到牛顿求解的精度：页顶默认参数下 $|q_M|=7.1\times10^{-15}$ ），所以含 $q_M$ 的项可以忽略。页顶默认参数下 $(g_h)_4=-\psi_{h,2}=-0.02009$ ， $(g_h)_5=-\psi_{h,3M-2}=-8.822\times10^{-4}$ 。图 14 把每一项画出来。

This is the heated case (the outlet constraints $u-0.35$ , $T-D_5$ of the manufactured case contain neither $\mu_{\rm ref}$ nor $Pr$ , the two terms with $\tau_M,q_M$ drop out, and $D_5$ comes from the temperature row, $\psi_{h,3M-1}$ ). $\overline R^{\rm raw}=M_{\rm bc}^T\psi_h$ has already zeroed the five constraint rows, so the two end faces do not appear; $q_M$ vanishes at the converged solution (to the accuracy of the Newton solve: $|q_M|=7.1\times10^{-15}$ at the page defaults), so the terms with $q_M$ are negligible. At the page defaults $(g_h)_4=-\psi_{h,2}=-0.02009$ and $(g_h)_5=-\psi_{h,3M-2}=-8.822\times10^{-4}$ . Figure 14 draws every term.

图 14 · 总梯度的每个分量从哪里来：第 9 节的反向扫描（种子

\overline R_h=\psi_h

、

\overline J_h=1

）按来源拆开（页顶默认参数，格点页，温度跟踪，

j_D=0

）。竖条：每个面对

\overline D_k

的贡献

\overline H_f^T\,\partial H_f/\partial D_k

（紫色，画在面位置）与每个控制体源项的贡献

Q_{ik}\overline S_{i,3}^h

（橙色，画在节点）；红色：被替换的五行约束的贡献（画在端点）。两个热源幅度只从源项得到；

\mu_{\rm ref}

与

Pr

从每个本构黏性通量得到；

D_4

与

D_5

各只从一处得到：入口温度约束行与出口应力约束行。每格的和等于

(g_h)_k

。

Figure 14 · Where each component of the total gradient comes from: the reverse sweep of Section 9 (seeds

\overline R_h=\psi_h

\overline J_h=1

) split by origin (page defaults, node page, temperature tracking,

j_D=0

). Bars: the contribution

\overline H_f^T\,\partial H_f/\partial D_k

of each face to

\overline D_k

(violet, drawn at the face) and the contribution

Q_{ik}\overline S_{i,3}^h

of each control volume's source (orange, drawn at the node); red: the contribution of the five replaced constraint rows (drawn at the ends). The two heat amplitudes come only from the sources;

\mu_{\rm ref}

and

Pr

from every constitutive viscous flux;

D_4

and

D_5

each from one place only: the inlet-temperature constraint row and the outlet-stress constraint row. Each panel sums to

(g_h)_k

下表把总梯度拆为显式目标项和残差反传项，并比较所选方向 $\dot D=e_k$ 的 tangent 结果 $\dot J_h$ 。若制造解源项在每次参数扰动后重新生成，比较的就不再是同一个函数；因此本页所有差分、复数计算和导数都使用基准模型冻结的制造源项。

The table separates direct objective and residual-reverse contributions and compares the tangent result $\dot J_h$ for the selected direction $\dot D=e_k$ . Regenerating manufactured forcing after every parameter perturbation would change the function being differentiated. All differences, complex calculations and derivatives therefore share the baseline model's frozen manufactured source.

第 13 节Section 13

有限差分与 Taylor 检验：重新求解，检查总导数

Finite differences and Taylor tests: re-solve and check total derivatives

有限差分必须重新解原始问题

Finite differences must re-solve the primal

为避免与空间网格尺度 $h$ 混淆，本节把参数差分步长记为 $\varepsilon$ 。选参数 $D_k$ 后分别形成 $D^\pm=D\pm\varepsilon e_k$ ，从基态作初值重新收敛原始方程，再计算 $\widehat J_h(D^\pm)$ 。只改变目标输入而固定流场，会验证偏导而不是总导数。

Use $\varepsilon$ for the parameter step to distinguish it from mesh spacing $h$ . For parameter $D_k$ , form $D^\pm=D\pm\varepsilon e_k$ , use the base state as an initial guess, reconverge the primal, and evaluate $\widehat J_h(D^\pm)$ . Perturbing only objective inputs while fixing the flow would test a partial derivative, not the total derivative.

\mathrm{FD}_k^+=\frac{\widehat J_h(D+\varepsilon e_k)-\widehat J_h(D)}\varepsilon,\quad \mathrm{FD}_k^-=\frac{\widehat J_h(D)-\widehat J_h(D-\varepsilon e_k)}\varepsilon,\quad \mathrm{FD}_k^0=\frac{\widehat J_h(D+\varepsilon e_k)-\widehat J_h(D-\varepsilon e_k)}{2\varepsilon}.

光滑且收敛充分时，前后差分截断误差通常为一阶，中心差分为二阶。步长过大受到非线性截断误差影响，过小则放大原始求解误差和浮点相消。这里参数有不同物理含义，统一绝对步长只是教学扫描方式，不等于所有参数的最佳步长。

For sufficiently converged smooth problems, forward and backward truncation errors are usually first order and centred errors second order. Large steps incur nonlinear truncation error; small steps amplify primal solve error and cancellation. A common absolute step is a teaching sweep, not necessarily optimal for parameters with different meanings.

E_1(\varepsilon)=|\widehat J_h(D+\varepsilon e_k)-\widehat J_h(D)|,\quad E_2(\varepsilon)=|\widehat J_h(D+\varepsilon e_k)-\widehat J_h(D)-\varepsilon(g_h)_k|.

一般非零一阶导数下， $E_1$ 为一阶、 $E_2$ 为二阶；零梯度目标的首项可能不同。下方可选所有参数、单个步长或完整步长扫描，并显示两个扰动问题的最终原始残差。判断时应寻找一段合理的收敛区间，而不是只取最小步长的一个数。

For a generic nonzero first derivative, $E_1$ is first order and $E_2$ second order; zero-gradient objectives can have different leading terms. Select any parameter, a single step or a complete sweep, with final primal residuals for both perturbed solves. Look for a credible convergence range rather than trusting one value at the smallest step.

四种误差要分开：空间离散误差由网格加密检查；非线性求解误差由原始残差检查；线性求解误差由 tangent/adjoint 真残差检查；求导实现由转置、差分与复数检查。点积相等只证明两条传播互为转置；共同遗漏同一项也可能通过（图 8）。

Separate four errors: spatial error through refinement; nonlinear solve error through primal residuals; linear solve error through true tangent/adjoint residuals; derivative implementation through transpose, finite differences and complex checks. A dot test alone only establishes paired transposes; jointly omitting a term may still pass (Figure 8).

第 14 节Section 14

复数步长：用虚部检验同一个离散函数

Complex step: check the same discrete function through its imaginary part

复数步长验证同一个离散函数

Complex step checks the same discrete function

(g_h)_k\approx\frac{\operatorname{Im}J_h(\mathbf U_h(D+\mathrm i\varepsilon e_k),D+\mathrm i\varepsilon e_k)}\varepsilon,\qquad\varepsilon=10^{-25}.

在当前正密度、正温度、无 limiter 的光滑分支，除法、平方根与 Sutherland 幂函数可沿邻近复数分支解析延拓。复数步长避免两个几乎相等的实目标相减，但不自动消除原始方程求解误差。复数状态必须满足复数残差，而不仅是把实解加上某个标签。

On the current positive-density, positive-temperature, unlimited smooth branch, division, square roots and the Sutherland power admit a nearby analytic complex continuation. Complex step avoids subtracting nearly equal real objectives, but does not automatically remove primal solve error. The complex state must satisfy the complex residual; attaching an imaginary part to the real solution is not enough.

本实现从已收敛实解出发，用独立复数算术评估整个残差与目标，再迭代校正缩放虚部残差。校正使用基态实 Jacobian 求解，这正是小虚部的一阶修正；最终目标仍由复数残差路径独立计算。它能检查该路径与 AD 导数是否一致，但不是对任意非解析代码的通用证明。

The implementation starts from the converged real state, evaluates the entire residual and objective using independent complex arithmetic, and corrects the scaled imaginary residual iteratively. Corrections use the real base Jacobian, appropriate to the small-imaginary-part linearization. The final objective is evaluated independently through the complex path. This checks consistency with AD derivatives, not arbitrary nonanalytic code.

r_I=\frac{\operatorname{Im}R_h(\mathbf U_h^{\mathbb C},D^{\mathbb C})}\varepsilon,\qquad A_h\Delta v=-r_I,\qquad\mathbf U_h^{\mathbb C}\leftarrow\mathbf U_h^{\mathbb C}+\mathrm i\varepsilon\Delta v.

复数路径不能对复数取绝对值、大小比较或裁剪后还声称保持解析性。当前有效性判断在实基态进行；如果以后加入 limiter 或激波捕捉，应单独分析分支，不直接沿用本节结论。运行检查后生成的表同时给出缩放虚部残差，防止很小的步长掩盖未收敛。

Complex absolute values, ordering comparisons or clipping generally destroy analyticity. Current admissibility checks use the real base state. Limiters or shock capture would require separate branch analysis rather than inheriting this result. The table produced by the check includes the scaled imaginary residual so that a tiny step cannot conceal nonconvergence.

第 15 节Section 15

精度阶：制造解上的一阶与二阶加密

Order of accuracy: first- and second-order refinement on a manufactured solution

本节实验：离散状态达到二阶，不等于目标或梯度也达到二阶。先运行制造解加密，读守恒量误差的实测阶，再到第 16 节比较梯度。

Experiment: Second-order states do not imply second-order objectives or gradients. Run the manufactured-solution refinement for the measured state rates, then compare gradients in Section 16.

空间精度与导数正确性是两件事

Spatial accuracy and derivative correctness are different tests

有限差分和点积检验在固定网格上检查导数；网格加密检查数值解逼近参考问题的速度。精确的一阶离散伴随仍然是一阶模型的导数，不会因为转置计算精确就自动变成二阶。

Finite differences and dot products test derivatives at fixed mesh; refinement tests how rapidly the numerical solution approaches a reference problem. An exact discrete adjoint of a first-order model remains that model's derivative; exact transposition does not confer second-order spatial accuracy.

u_*(x)=0.35+0.03\sin(\pi x)+0.01\sin(2\pi x),\quad T_*(x)=1+0.05\sin(\pi x)+0.02\sin(2\pi x),\quad\rho_*(x)=0.35/u_*(x).

把这些解析函数代入物理总通量并作解析空间求导，构造独立的制造源项；质量通量恒为 0.35，质量源为零。这一步是生成参考问题，不是用连续伴随代替离散求导。数值状态仍由本页残差求解。参考源项包含物性随温度的变化，不能把变物性情形当作常系数二阶导数。

Substitute these analytic functions into the physical net flux and differentiate analytically in space to manufacture an independent source. Mass flux is constant at 0.35, giving zero mass source. This constructs the reference problem; it does not substitute continuous adjoints for discrete differentiation. The numerical state is still obtained from this page's residual. Variable transport derivatives belong in the manufactured forcing and cannot be replaced by constant-coefficient second derivatives.

E_h=\left(\sum_iV_i\|U_i-U_{*,i}\|_2^2\right)^{1/2},\qquad p_h=\frac{\log(E_h/E_{h/2})}{\log2}.

格点参考 $U_{*,i}$ 是解析节点值。端点误差按半控制体体积计权。

The node reference $U_{*,i}$ is the analytic point value. Endpoint errors use half-volume weights.

下方独立运行制造解加密，分别计算一阶、二阶误差和实测阶数，保留当前物性选择。有限网格上的斜率受对流、扩散和边界误差共同影响；参考斜率线只是视觉标尺。应确认原始残差远小于空间误差，并观察连续多个加密层级，不能把单个粗网格斜率当作阶数证明。

The panel independently refines the manufactured case for first and second order, retaining current transport settings. Finite-grid rates reflect convection, diffusion and boundary errors together; reference slopes are visual guides. Check that primal residuals are much smaller than spatial errors and inspect several refinement levels rather than treating one coarse-grid slope as proof.

固定基准补充： $\mu_{\mathrm{ref}}$ =0.15、Pr=0.72、常物性；下表是离线实际计算的各级网格误差与相邻加密斜率，不随顶部参数改变。这些网格还没有完全进入渐近区：一阶斜率从 0.4–0.5 逐渐升向 1，二阶斜率高于 2 并随加密回落。它们只说明这个网格范围内误差的下降速度，不能单独当作渐近阶的证明。

Fixed supplementary baseline: $\mu_{\mathrm{ref}}$ =0.15, Pr=0.72, constant transport. The table gives measured offline errors on every grid and the rates between neighbouring grids; it does not follow the top controls. These grids are not yet fully asymptotic: first-order slopes rise from 0.4–0.5 towards 1, and second-order slopes exceed 2 and fall back under refinement. They describe how fast the error drops over this range, not a proof of the asymptotic order.

阶数order	N	误差error	实测阶measured rate
1	8	1.16842e-1	—
1	16	8.89070e-2	0.3942
1	32	6.06293e-2	0.5523
1	64	3.79622e-2	0.6755
1	128	2.21712e-2	0.7759
2	8	2.07603e-2	—
2	16	4.39520e-3	2.2398
2	32	9.52454e-4	2.2062
2	64	2.16414e-4	2.1379
2	128	5.11940e-5	2.0797

第 16 节Section 16

伴随一致性：离散梯度与连续参考分开检查

Adjoint consistency: discrete gradients and the continuous reference checked separately

先说明本节比较的对象

Define what is being compared

离散点积测试验证实现是否互为转置；差分和复数步长验证是否对当前离散目标求了完整导数。本节问的是另一个问题：网格加密时，目标与总设计梯度是否逼近同一个物理边值问题。连续参考的完整推导和独立实现放在附录 A，正文不使用连续伴随微分方程求解 $\psi_h$ 。

Discrete dot tests check transposition; finite differences and complex step check the complete derivative of the current discrete objective. This section asks a different question: do objective and total design gradients approach the same physical boundary-value problem under refinement? Appendix A contains the continuous derivation and independent implementation. The main discrete solve never uses a continuous adjoint differential equation to obtain $\psi_h$ .

连续参考先利用恒定质量流率消去密度，采用速度、温度、总动量通量和总能量通量四个变量。其四个连续乘子记作 $\lambda^{\rm c}(x)$ ；正文的 $\psi_h$ 有 $3M$ 个离散行乘子。两者的变量、维数和边界处理不同；本节下方的对照图比较目标与同一参数的总导数。乘子本身只有动量行与能量行可以逐点对应，推导如下。

The continuous reference eliminates density using constant mass flux and uses velocity, temperature, total momentum flux and total energy flux as four variables. Its four continuous multipliers are denoted $\lambda^{\rm c}(x)$ , while $\psi_h$ contains $3M$ discrete row multipliers. Variables, dimensions and boundary treatments differ; the comparison panel below compares objectives and total derivatives for the same parameter. Of the multipliers themselves, only the momentum and energy rows correspond pointwise, as derived next.

离散乘子与连续乘子的对应

Matching the discrete and continuous multipliers

离散残差是控制体上的积分平衡。动量行和能量行的数值总通量，正是附录 A.2 连续参考所用的两个变量：总动量通量 $\Pi=m_*u+p-\tau$ 与总能量通量 $\mathcal H=m_*(c_pT+u^2/2)-u\tau+q$ 。所以这两行就是 $\Pi'=S_m$ 、 $\mathcal H'=S_E$ 在控制体上的积分：

The discrete residual is an integrated balance over a control volume. The numerical total fluxes of the momentum and energy rows are exactly the two variables of the continuous reference in Appendix A.2: the total momentum flux $\Pi=m_*u+p-\tau$ and the total energy flux $\mathcal H=m_*(c_pT+u^2/2)-u\tau+q$ . These two rows are therefore $\Pi'=S_m$ and $\mathcal H'=S_E$ integrated over a control volume:

\begin{gathered}R_{i,2}=\Pi_{i+1/2}-\Pi_{i-1/2}-\int_{V_i}S_m\,dx\approx\int_{V_i}\bigl(\Pi'-S_m\bigr)dx,\\R_{i,3}=\mathcal H_{i+1/2}-\mathcal H_{i-1/2}-\int_{V_i}S_E\,dx\approx\int_{V_i}\bigl(\mathcal H'-S_E\bigr)dx.\end{gathered}

比较两个 Lagrangian 中的对应项。离散约定是 $+\psi^TR$ ，附录 A.3 的连续约定是 $\lambda^T(f-z')$ ，而 $f$ 的 $\Pi$ 、 $\mathcal H$ 分量正是 $S_m$ 、 $S_E$ ：

Compare the corresponding terms of the two Lagrangians. The discrete convention is $+\psi^TR$ and the continuous one in Appendix A.3 is $\lambda^T(f-z')$ , whose $\Pi$ and $\mathcal H$ components of $f$ are exactly $S_m$ and $S_E$ :

\begin{gathered}\sum_i\psi_{i,2}R_{i,2}\approx\sum_i\psi_{i,2}\int_{V_i}\bigl(\Pi'-S_m\bigr)dx,\\\int_0^1\lambda_\Pi\bigl(S_m-\Pi'\bigr)dx=-\int_0^1\lambda_\Pi\bigl(\Pi'-S_m\bigr)dx.\end{gathered}

同一个约束的乘子只差这个符号，因此在光滑区域 $\psi_{i,2}\approx-\lambda_\Pi(x_i)$ 、 $\psi_{i,3}\approx-\lambda_{\mathcal H}(x_i)$ 。残差没有除以控制体体积，所以两组乘子量级相同。另外两个连续乘子没有离散对应： $\lambda_u$ 、 $\lambda_T$ 属于本构关系 $u'=\tau/(4\mu/3)$ 与 $T'$ 的热流关系，有限体积在黏性通量内部计算它们，不是残差行；质量行也没有对应，因为连续参考用固定的 $m$ 消去了密度，没有质量方程。

Multipliers of the same constraint differ only by this sign, so in smooth regions $\psi_{i,2}\approx-\lambda_\Pi(x_i)$ and $\psi_{i,3}\approx-\lambda_{\mathcal H}(x_i)$ . The residual is not divided by the control volume, so both sets have the same magnitude. The other two continuous multipliers have no discrete counterpart: $\lambda_u$ and $\lambda_T$ belong to the constitutive relations $u'=\tau/(4\mu/3)$ and the heat-flux relation for $T'$ , which the finite-volume scheme evaluates inside the viscous flux rather than as residual rows; the mass row has no counterpart either, because the continuous reference eliminates density with the fixed $m$ and has no mass equation.

图 15 · 本页格点格式在 N = 24 时的动量、能量行乘子（点）与附录 A 的连续乘子取负号（曲线），目标为默认的温度跟踪。一阶点整体偏离曲线，是 O(h) 的内部误差。格点入口三行与出口两行是约束，其乘子不画；紧邻出口约束的节点明显偏离曲线。

Figure 15 · Momentum and energy row multipliers (points) of this page's node scheme at N = 24, against the negated continuous multipliers of Appendix A (curves), for the default temperature-tracking objective. The first-order points sit off the curves throughout: that is their O(h) interior error. The three inlet rows and two outlet rows of the node layout are constraints and are not drawn; the node next to the outlet constraints leaves the curve visibly.

下表对本页格式逐级加密，比较 $\max|\psi_{i,2}+\lambda_\Pi(x_i)|$ 与 $\max|\psi_{i,3}+\lambda_{\mathcal H}(x_i)|$ ；“内部”取 $0.15<x<0.85$ ，“全域”取全部场行及其位置。连续乘子由附录 A.3–A.4 的独立 ODE 伴随在 768 步 RK4 上求得，每个网格点都是它的节点。

The table refines this page's scheme and compares $\max|\psi_{i,2}+\lambda_\Pi(x_i)|$ and $\max|\psi_{i,3}+\lambda_{\mathcal H}(x_i)|$ ; "interior" means $0.15<x<0.85$ , "global" all field rows, with the location of the maximum. The continuous multipliers come from the independent ODE adjoint of Appendices A.3–A.4 on 768 RK4 steps, and every mesh point is one of its nodes.

阶数Order	N	内部动量Interior momentum	阶Rate	内部能量Interior energy	阶Rate	全域动量Global momentum	全域能量Global energy
二阶second order	12	3.07e−4	—	6.85e−5	—	1.64e−3 (x = 0.917)	1.68e−3 (x = 0.917)
	24	8.88e−5	1.79	2.18e−5	1.65	1.35e−3 (x = 0.958)	1.90e−3 (x = 0.958)
	48	2.80e−5	1.67	6.25e−6	1.80	1.21e−3 (x = 0.979)	2.03e−3 (x = 0.979)
	96	7.93e−6	1.82	1.72e−6	1.87	1.15e−3 (x = 0.990)	2.10e−3 (x = 0.990)
一阶first order	12	2.03e−3	—	8.72e−4	—	3.68e−3 (x = 0.917)	2.51e−3 (x = 0.917)
	24	1.43e−3	0.51	4.68e−4	0.90	2.82e−3 (x = 0.958)	2.52e−3 (x = 0.958)
	48	8.96e−4	0.67	2.43e−4	0.95	2.12e−3 (x = 0.979)	2.41e−3 (x = 0.979)
	96	5.19e−4	0.79	1.23e−4	0.99	1.67e−3 (x = 0.990)	2.32e−3 (x = 0.990)

内部按格式的阶收敛（二阶约 2，一阶逐渐接近 1）。全域最大偏差都出现在紧邻出口约束的节点 $x=1-h$ ，而且不随加密减小：二阶从 $N=12$ 到 $96$ ，动量只从 $1.6\times10^{-3}$ 降到 $1.2\times10^{-3}$ ，能量反而从 $1.7\times10^{-3}$ 增到 $2.1\times10^{-3}$ ，与乘子本身同一量级。出口的动量行和能量行已换成约束，相邻节点的平衡行承担了边界的作用，它的乘子不是连续场的取值。这与一维 Burgers 格点页入口的伴随边界层类似：总导数仍然收敛，六个参数导数的最大相对误差在二阶格点上从 $3.2\times10^{-1}$ 降到 $1.0\times10^{-2}$ （最后一级实测阶 1.81），但逐点比较乘子时必须排除这些边界点。

The interior converges at the order of the scheme (about 2 at second order, approaching 1 at first order). The global maximum always sits at the node next to the outlet constraints, $x=1-h$ , and does not decrease under refinement: at second order from $N=12$ to $96$ momentum falls only from $1.6\times10^{-3}$ to $1.2\times10^{-3}$ and energy even grows from $1.7\times10^{-3}$ to $2.1\times10^{-3}$ , the same magnitude as the multipliers themselves. The outlet momentum and energy rows have been replaced by constraints, the balance rows of the neighbouring node take over the boundary role, and its multipliers are not values of the continuous field. This resembles the inlet adjoint layer of the 1-D Burgers node page: total derivatives still converge, the maximum relative error of the six parameter derivatives falling from $3.2\times10^{-1}$ to $1.0\times10^{-2}$ at second order (last measured rate 1.81), but a pointwise comparison of multipliers must exclude these boundary points.

残差缩放会改变离散乘子的数值

Residual scaling changes discrete multiplier values

\widetilde R_h=\Gamma^{-1}R_h,\quad\widetilde A_h=\Gamma^{-1}A_h,\quad\widetilde\psi_h=\Gamma^T\psi_h,\quad\widetilde A_h^T\widetilde\psi_h=-j_U.

以上在精确收敛状态成立， $\Gamma$ 是可逆行尺度矩阵。把平衡行除以控制体体积时，其乘子相应乘以体积；总设计梯度不变。即使 $\Gamma$ 随参数改变，新增项也乘以零残差。边界约束行需要按自己的尺度处理，不能不加区分地套用体积权重。

These identities hold at exact convergence for an invertible row-scaling matrix $\Gamma$ . Dividing balance rows by volume multiplies their corresponding multipliers by volume, preserving the total design gradient. Parameter-dependent scaling adds terms multiplied by the zero residual. Boundary constraint rows need their own scaling; volume weights cannot be applied indiscriminately.

本页转置使用普通数组点积，不是带体积权的函数空间伴随。如果引入权矩阵，伴随算子需相应变换。连续参考也有 ODE 积分与射击误差；先提高参考分辨率，再比较离散加密误差，避免把参考误差平台误读为空间离散不一致。

This page uses ordinary array dot products, not a volume-weighted function-space adjoint. Introducing weight matrices requires transforming the adjoint operator. The continuous reference also has ODE integration and shooting errors; refine it before interpreting a discrete-refinement plateau as spatial inconsistency.

粗网格梯度误差可以不单调

Gradient error need not decrease monotonically on coarse grids

本节的表是格心基准，两页相同；格点页的梯度可在附录 A 按当前设置与连续参考比较。

The table in this section is the cell baseline, identical on both pages; on the node page, compare gradients with the continuous reference for the current settings in Appendix A.

固定基准为格心、一阶重构、局部加热 $a_0=0.03$ 、常物性 $\mu_{\rm ref}=0.15,Pr=0.72$ 、温度跟踪目标；比较对出口总应力参数 $D_5$ 的总梯度，连续参考使用 512 个 ODE 步。离线实测绝对误差如下，不随顶部控件改变：

For a fixed baseline, use cells, first-order reconstruction, local heating $a_0=0.03$ , constant transport $\mu_{\rm ref}=0.15,Pr=0.72$ , and temperature tracking. Compare the total derivative with respect to outlet total stress $D_5$ against a 512-step ODE reference. These measured offline absolute errors do not follow the top controls:

N	$\|(g_h)_5-(g^{\rm c})_5\|$
8	1.57965e−4
16	2.55325e−4
32	2.31363e−4
64	1.64320e−4
128	1.01779e−4

误差先增大，再随加密下降。该数据说明，粗网格上的偶然误差抵消不能当作渐近精度；它本身也不能证明伴随发散或连续不一致。第 13 节检查固定网格上的求导，第 15 节检查状态精度，本节检查目标与梯度的连续对照，三者不可互相替代。在格心页设置上述基准、选择出口参数并运行加密，可重现前四行；128 网格数据是单独的离线补充。

Error first increases, then falls with refinement. Accidental coarse-grid error cancellation does not establish asymptotic accuracy; this sequence alone proves neither adjoint divergence nor continuum inconsistency. Section 13 checks differentiation at a fixed mesh, Section 15 checks state accuracy, and this section compares objectives and gradients with a continuum reference. These tests are not interchangeable. On the cell page, select the stated baseline and outlet parameter, then run refinement to reproduce the first four rows; N=128 is a separate offline supplement.

第 17 节Section 17

GS 固定点与 GS 预处理 GMRES

GS fixed point and GS-preconditioned GMRES

GS 缺陷校正解的是精确系统

GS defect correction targets the exact system

设目标线性系统为 $Kx=b$ ：tangent 时 $K=A_h$ ，adjoint 时 $K=A_h^T$ 。每轮先用精确算子求残差，再用较便宜的 GS 映射修正。不能只把 $P_hx=b$ 解得很准，就宣称解好了原系统。

Write the target system as $Kx=b$ : use $K=A_h$ for tangents and $K=A_h^T$ for adjoints. Each iteration computes the residual with the exact operator, then applies the cheaper GS map. Accurately solving $P_hx=b$ alone does not solve the original system.

r^{(k)}=b-Kx^{(k)},\qquad x^{(k+1)}=x^{(k)}+\omega B_Kr^{(k)},\quad\omega=0.65,\qquad B_K=\begin{cases}B_h&K=A_h,\\B_h^T&K=A_h^T.\end{cases}

e^{(k+1)}=(I-\omega B_KK)e^{(k)},\qquad\rho(I-\omega B_KK)<1\quad\Longrightarrow\quad e^{(k)}\to0.

这是固定线性系统下的收敛条件；当前有限参数试验不是一般收敛定理。黏性耦合、重构和边界可能使一次 GS 近似不足。预算用尽表示本次未收敛，页面仍报告精确算子的真残差。

This is a convergence condition for the fixed linear system, not a general theorem inferred from the exposed parameter cases. Viscous coupling, reconstruction and boundaries can make one GS approximation insufficient. Exhausting the budget means this run has not converged; the page still reports the exact-system residual.

转置 GS 必须反向扫描

Transposed GS requires a reverse sweep

\begin{gathered}B_h^T=(\operatorname{diag}(P_h)^T+L_P^T)^{-1},\\(B_h^Tr)_i=P_{ii}^{-T}\left(r_i-\sum_{j>i}P_{ji}^T(B_h^Tr)_j\right),\quad i=M-1,\ldots,0.\end{gathered}

注意求和用 $P_{ji}^T$ ，不是 $P_{ij}^T$ 。先处理后面的块，才能读到已经求出的分量；每个对角块也要转置。直接对原 GS 程序改循环方向但不转置块或累加路径，不会得到正确的伴随预处理。

The sum uses $P_{ji}^T$ , not $P_{ij}^T$ . Later blocks must be solved first so their components are available, and every diagonal block must be transposed. Reversing only the loop direction without transposing blocks and accumulation paths does not produce the adjoint preconditioner.

GMRES：用 GS 作右预处理

GMRES with GS as a right preconditioner

\beta=\|b\|_2,\quad KB_Ky=b,\quad x=B_Ky,\qquad KB_KV_m=V_{m+1}\underline H_m,\quad y_m=\arg\min_y\|\beta e_1-\underline H_my\|_2,\quad x_m=B_KV_my_m.

GMRES 不再只沿单个残差修正，而是在多个 Krylov 方向张成的空间中最小化残差。 $\underline H_m$ 是 Arnoldi 过程的 $(m+1)\times m$ 上 Hessenberg 矩阵（上横线在本页留给反向模式导数，所以这里用下横线）。本实现每个新方向先做一次固定的 GS 预处理，再调用精确乘积；保存预处理后的方向来恢复解。Arnoldi 正交化和小型最小二乘代价也属于求解成本。

GMRES minimizes residuals over a space of multiple Krylov directions instead of taking only one residual correction at a time. Here $\underline H_m$ is the $(m+1)\times m$ upper Hessenberg matrix of the Arnoldi process (underlined, because this page reserves the overbar for reverse-mode derivatives). Each new direction receives one fixed GS preconditioning application followed by an exact operator product. Preconditioned directions are stored to recover the solution. Arnoldi orthogonalization and the small least-squares problem also contribute to cost.

b = -jU; x = 0
# GS correction
repeat:
    r = b - exact_AT(x)
    stop if norm(r) < tolerance
    x += 0.65 * reverse_block_GS(r)
# GMRES
for each Arnoldi direction v:
    z = reverse_block_GS(v)
    w = exact_AT(z)
    orthogonalize w; update least-squares residual
    recover x from stored z directions
    confirm norm(b - exact_AT(x))

下方在同一收敛流场上比较 GS、GMRES 和 LU，计时包含各自准备。GMRES 有基于正交化的残差估计，也计算恢复解的真残差；最终判据应看真残差。小网格上 LU 可能更快，附录 C 明确展示其显式组装与分解。这里不把线性迭代次数与外层 Newton 次数混为一谈。

The panel compares GS, GMRES and LU on the same converged flow, including each method's preparation. GMRES has an orthogonalization-based residual estimate and also evaluates the recovered solution's true residual; final acceptance should use the latter. LU can be faster on this small mesh; Appendix C explicitly shows assembly and factorization. Linear iterations are distinct from outer Newton iterations.

第 18 节Section 18

总结：从守恒状态到经过验证的设计梯度

Summary: from conservative states to verified design gradients

把整条计算链连起来

Connect the entire computational chain

本页从有限个守恒状态出发：重构决定对流面状态，原始量转换和差分决定黏性面数据，物性决定应力和热流，面装配与边界覆盖定义最终残差。目标是另一个有限维标量函数。任何一条状态或参数依赖漏掉，最终伴随就会求解另一套导数问题。

The lesson starts with finitely many conservative states. Reconstruction defines convective traces; primitive conversion and differences define viscous face data; transport defines stress and heat flux; face assembly and boundary overwrites define the final residual. The objective is another finite-dimensional scalar function. Omitting any state or parameter dependency changes the derivative problem solved by the adjoint.

R_h(\mathbf U_h,D)=0,\qquad A_hv=-C_h\dot D,\qquad A_h^T\psi_h=-j_U,\qquad g_h=j_D+C_h^T\psi_h.

这四个式子分别是原始方程、设计 tangent（ $v\equiv\dot{\mathbf U}_h$ 沿设计方向 $\dot D$ ）、离散 adjoint 和总设计梯度（以 $\overline R_h=\psi_h$ 、 $\overline J_h=1$ 反向扫描得到的 $\overline D$ ）。matrix-free 只改变乘积的实现方式；GS、GMRES、LU 只改变线性系统的求解方式；空间一阶／二阶则会改变残差及其精确 Jacobian。三者不能互相替代。

These four equations are the primal, the design tangent ( $v\equiv\dot{\mathbf U}_h$ along the design direction $\dot D$ ), the discrete adjoint and the total design gradient (the $\overline D$ of the reverse sweep seeded with $\overline R_h=\psi_h$ , $\overline J_h=1$ ). Matrix-free evaluation changes how products are computed. GS, GMRES and LU change how a linear system is solved. First- versus second-order space discretization changes the residual and its exact Jacobian. These choices are not interchangeable.

建议按层核对：先看局部面求导与点积，再看全局边界掩码与转置，再看线性真残差，然后做完整重算的差分和复数验证，最后做空间加密与独立连续参考对照。每个检查回答不同问题，没有一个单独检查能替代全部验证。

Check the layers in order: local face derivatives and dot products, global boundary masking and transposition, true linear residuals, fully re-solved finite differences and complex step, then spatial refinement and independent continuous comparison. Each answers a different question; none replaces all the others.

当前实现是固定一维网格上的光滑层流基础算例。走向二维非结构网格时，需要扩展多维梯度重构、面法向、黏性边界和几何求导；走向 RANS 时还需湍流模型及其完整耦合导数。这些是后续目标，不应与本页已经验证的范围混淆。

The present implementation is a smooth laminar foundation on a fixed one-dimensional mesh. Two-dimensional unstructured grids require multidimensional gradients, face normals, viscous boundaries and geometry derivatives. RANS additionally requires turbulence equations and their complete coupled derivatives. These are subsequent goals, distinct from the scope verified here.

附录 A：独立连续参考与完整伴随

Appendix A: independent continuous reference and complete adjoint

进阶教学 · 选读Advanced teaching · optional

这里是已有连续伴随理论在一维光滑层流上的推导与独立验证，不因推导较长而归为研究成果。质量方程没有扩散，黏性系数矩阵因而奇异，所以改用四维约化 ODE。出现分辨不足的边界层或改变边界类型时，需要重新检查推导的前提；不能仅由离散／连续不一致判定连续问题不适定。

This is established continuous-adjoint theory derived and independently checked for smooth one-dimensional laminar flow, not a research result because the derivation is lengthy. The mass equation has no diffusion, so the viscous coefficient matrix is singular and a reduced four-variable ODE is used. Under-resolved layers or different boundary types require revisiting its assumptions; discrete/continuous disagreement alone does not establish ill-posedness.

A.1 本附录才使用连续场和空间微分

A.1 Continuous fields and spatial derivatives belong here

本附录构造独立连续参考。用 $U^{\rm c}(x)$ 表示连续守恒场，用 $\chi^{\rm c}(x)$ 表示三分量形式连续乘子；约化 ODE 的四分量乘子另记 $\lambda^{\rm c}(x)$ 。正文的 $\mathbf U_h,R_h,A_h,\psi_h$ 仍专指有限维离散量。本附录中的积分、空间导数及分部积分不参与正文离散转置的定义。

This appendix constructs an independent continuous reference. Denote the conservative field by $U^{\rm c}(x)$ and its three-component formal continuous multiplier by $\chi^{\rm c}(x)$ . The four-component multiplier of the reduced ODE is $\lambda^{\rm c}(x)$ . The main text's $\mathbf U_h,R_h,A_h,\psi_h$ remain finite-dimensional discrete quantities. Integrals, spatial derivatives and integration by parts here do not define the main discrete transpose.

\mathcal R^{\rm c}=\frac{d}{dx}F(U^{\rm c},U_x^{\rm c},D)-S=0,\quad F=F^c-F^v,\quad \tau=\frac43\mu u_x,\quad q=-\kappa T_x.

先考虑体积分目标密度 $\ell$ 不依赖状态空间导数的情形；端点目标保留为 $\phi_0,\phi_1$ 。把残差以加号接入： $\mathcal L^\chi=\int_0^1\ell\,dx+\phi_0+\phi_1+\int_0^1(\chi^{\rm c})^T\mathcal R^{\rm c}\,dx$ ；用通量 Jacobian $F_U,F_{U_x}$ 线性化残差，并两次分部积分。以下梯度下标表示行偏导，连续空间导数明确写作 $d/dx$ 或下标 x。

First let the volume objective density $\ell$ be independent of spatial state derivatives, retaining endpoint objectives $\phi_0,\phi_1$ . Adjoin the residual with a plus sign, $\mathcal L^\chi=\int_0^1\ell\,dx+\phi_0+\phi_1+\int_0^1(\chi^{\rm c})^T\mathcal R^{\rm c}\,dx$ ; linearize the residual with the flux Jacobians $F_U,F_{U_x}$ and integrate by parts twice. Subscripts below (as in $F_U$ , $F_{U_x}$ , $F_D$ ) denote partial derivatives, written as row vectors; spatial derivatives are explicitly written as $d/dx$ or subscript x.

\delta\mathcal R^{\rm c}=\frac d{dx}(F_U\delta U^{\rm c}+F_{U_x}\delta U_x^{\rm c}+F_D\delta D)-S_U\delta U^{\rm c}-S_D\delta D,

\ell_U^T-F_U^T\chi_x^{\rm c}+\frac d{dx}(F_{U_x}^T\chi_x^{\rm c})-S_U^T\chi^{\rm c}=0,

\left[(\chi^{\rm c})^T(F_U\delta U^{\rm c}+F_{U_x}\delta U_x^{\rm c}+F_D\delta D)-(\chi_x^{\rm c})^TF_{U_x}\delta U^{\rm c}\right]_0^1+\delta\phi_0+\delta\phi_1.

这给出形式体方程与边界项，还没有自动给出可求解的边值问题。黏性通量的质量分量为零， $F_{U_x}$ 奇异，不能把它当作三个独立二阶方程的可逆最高阶系数。代码实际采用下面的约化系统，并对该系统推导连续伴随。

These are formal interior equations and boundary terms, not yet an automatically solvable boundary-value problem. The viscous mass flux is zero, making $F_{U_x}$ singular; it is not an invertible highest-order coefficient for three independent second-order equations. The implemented continuous reference instead uses the reduced system below and its own derived adjoint.

参数变化还产生体积分 $\int_0^1[\ell_D^T-F_D^T\chi_x^{\rm c}-S_D^T\chi^{\rm c}]\,dx$ ；端点还含上式的 $F_D\delta D$ 与边界约束变化。尤其黏度和 Prandtl 数影响通量，不能只对源项求导。下面约化系统的 $f_D$ 和端点约束导数统一包含这些依赖。

Parameter variation also gives the volume contribution $\int_0^1[\ell_D^T-F_D^T\chi_x^{\rm c}-S_D^T\chi^{\rm c}]\,dx$ ; endpoints also contain the displayed $F_D\delta D$ and boundary-constraint variations. Viscosity and Prandtl number affect fluxes, so differentiating only the source is insufficient. The reduced system below includes these dependencies through $f_D$ and endpoint derivatives.

逐步推导： $F_{U_x}$ 的秩与约化后的阶数

Step by step: the rank of $F_{U_x}$ and the order of the reduced problem

黏性通量是 $F^v=(0,\tau,u\tau-q)^T$ ， $\tau=\frac43\mu u_x$ ， $q=-\kappa T_x$ 。梯度对 $U_x$ 是线性的： $u_x=u_UU_x$ ， $T_x=T_UU_x$ ，其中 $U=(\rho,m,e)^T$ 、 $E=e/\rho$ ：

The viscous flux is $F^v=(0,\tau,u\tau-q)^T$ with $\tau=\frac43\mu u_x$ and $q=-\kappa T_x$ . The gradients are linear in $U_x$ : $u_x=u_UU_x$ and $T_x=T_UU_x$ , where $U=(\rho,m,e)^T$ and $E=e/\rho$ :

\begin{gathered}u_U=\Bigl(-\frac u\rho,\ \frac1\rho,\ 0\Bigr),\qquad T_U=\frac1{\rho c_v}\bigl(u^2-E,\ -u,\ 1\bigr),\\ F_{U_x}=-F^v_{U_x}=-\begin{bmatrix}0\\ \frac43\mu\,u_U\\ \frac43\mu\,u\,u_U+\kappa\,T_U\end{bmatrix}.\end{gathered}

$F_{U_x}$ 的第一行为零；后两行中 $u_U$ 的第三个分量为零而 $T_U$ 的不为零，两者线性无关，所以 $\operatorname{rank}\,F_{U_x}=2$ 。于是三条守恒律的阶数是：质量方程只含一阶导数，动量、能量方程分别对 $u$ 、 $T$ 是二阶。总阶数正好等于原始问题的边界条件个数：

The first row of $F_{U_x}$ is zero; in the other two rows the third component of $u_U$ is zero while that of $T_U$ is not, so they are linearly independent and $\operatorname{rank}\,F_{U_x}=2$ . The mass equation therefore contains only first derivatives, while the momentum and energy equations are second order in $u$ and $T$ . The total order equals the number of primal boundary conditions:

\underbrace{1}_{(\rho u)'}+\underbrace{2}_{u''}+\underbrace{2}_{T''}=5=\underbrace{3}_{\rho(0),\,u(0),\,T(0)}+\underbrace{2}_{g_1=0}.

质量方程可以直接积分： $\rho u\equiv m_*$ ，积分常数用掉入口的一个条件（本页 $\rho(0)=1$ 、 $u(0)=m_*$ ）。代入 $\rho=m_*/u$ 后剩下 $u$ 、 $T$ 两个二阶方程；把 $u'$ 、 $T'$ 换成通量 $\Pi$ 、 $\mathcal H$ ，就是 A.2 的四维一阶系统，入口、出口各两个条件。格点页的五个约束行（入口三行、出口两行）正是这 $3+2$ 个条件。

The mass equation integrates directly to $\rho u\equiv m_*$ ; the integration constant uses one inlet condition (here $\rho(0)=1$ and $u(0)=m_*$ ). With $\rho=m_*/u$ , two second-order equations in $u$ and $T$ remain; replacing $u'$ and $T'$ by the fluxes $\Pi$ and $\mathcal H$ gives the four-dimensional first-order system of A.2, with two conditions at each end. The five constraint rows of the node page (three at the inlet, two at the outlet) are exactly these $3+2$ conditions.

A.2 消去恒定质量流率，建立四维 ODE

A.2 Eliminate constant mass flux to form a four-dimensional ODE

m_*=0.35,\quad\rho=m_*/u,\quad z=(u,T,\Pi,\mathcal H)^T,\quad \Pi=m_*u+p-\tau,\quad\mathcal H=m_*(c_pT+u^2/2)-u\tau+q.

\tau=m_*u+m_*T/u-\Pi,\quad q=\mathcal H-m_*(c_pT+u^2/2)+u\tau,\quad \frac{dz}{dx}=f(z,D,x)=\begin{bmatrix}\tau/(4\mu/3)\\-q/\kappa\\S_m\\S_E\end{bmatrix}.

$\Pi,\mathcal H$ 是总动量与总能量通量。入口给定 $u(0)=m_*,T(0)=D_4$ ，未知入口通量为两个射击变量。出口加热条件为 $g_1=(\Pi-m_*u-D_5,q)^T=0$ ；制造解条件为 $g_1=(u-m_*,T-D_5)^T=0$ 。本页不把入口密度或速度作为设计变量，因此质量流率固定。

The variables $\Pi,\mathcal H$ are total momentum and energy flux. Inlet conditions are $u(0)=m_*,T(0)=D_4$ , leaving two unknown inlet fluxes as shooting variables. The heated outlet imposes $g_1=(\Pi-m_*u-D_5,q)^T=0$ ; the manufactured outlet imposes $g_1=(u-m_*,T-D_5)^T=0$ . Inlet density and velocity are not design variables here, so mass flux stays fixed.

用 RK4 前向积分，再用两变量 Newton 修正入口通量以满足出口条件。射击 Jacobian 来自 ODE 对初值的变分方程。保存轨迹及右端值，使用三次 Hermite 插值供反向伴随积分取值。较强加热使用从零热源逐步增加的延拓；积分越出正值亚声速分支时报告失败，不裁剪状态。

Integrate forward with RK4 and use a two-variable Newton update of inlet fluxes to satisfy outlet conditions. Shooting derivatives come from ODE initial-value variations. Stored trajectories and right-hand sides provide cubic Hermite interpolation for backward adjoint integration. Stronger heating uses continuation from zero heating. Leaving the positive subsonic branch reports failure rather than clipping states.

A.3 连续目标、积分号与符号约定

A.3 Continuous objectives, integrals and signs

J^{\rm c}=\int_0^1\ell(z,D,x)\,dx+\phi_0(z(0),D)+\phi_1(z(1),D),

\mathcal L^{\rm c}=\int_0^1[\ell+(\lambda^{\rm c})^T(f-z_x)]\,dx+\phi_0+\phi_1+\eta_0^Tg_0+\eta_1^Tg_1,\quad g_0=(u(0)-m_*,T(0)-D_4)^T.

温度跟踪使用 $\ell=(T-T^*)^2/2$ 、零端点目标；出口温度使用 $\phi_1=T(1)$ ；入口热流使用 $\phi_0=q(0)$ 。这里乘的是 $f-z_x$ ，与正文的加号离散残差约定不是同一个对象。对 $-(\lambda^{\rm c})^T\delta z_x$ 作一次分部积分，分别收集体内、入口和出口的系数，得到：

Tracking uses $\ell=(T-T^*)^2/2$ with zero endpoint objectives; outlet temperature uses $\phi_1=T(1)$ and inlet heat flux uses $\phi_0=q(0)$ . The multiplier acts on $f-z_x$ , a different object from the main discrete residual convention. Integrating $-(\lambda^{\rm c})^T\delta z_x$ by parts once and collecting interior and endpoint coefficients gives:

\frac{d\lambda^{\rm c}}{dx}=-f_z^T\lambda^{\rm c}-\ell_z^T,

\lambda^{\rm c}(0)+\phi_{0,z}^T+g_{0,z}^T\eta_0=0,\qquad-\lambda^{\rm c}(1)+\phi_{1,z}^T+g_{1,z}^T\eta_1=0.

逐步推导：伴随在每一端给两个条件

Step by step: two adjoint conditions at each end

A.3 的两个端点方程各有四个分量，而 $\eta_0$ 、 $\eta_1$ 各只有两个未知数；消去它们，每端剩两个只含 $\lambda^{\rm c}$ 的条件。入口 $g_0=(u(0)-m_*,T(0)-D_4)^T$ ， $g_{0,z}$ 的两行是 $e_u^T$ 、 $e_T^T$ ，所以入口方程的后两个分量不含 $\eta_0$ ：

Each endpoint equation of A.3 has four components, while $\eta_0$ and $\eta_1$ have only two unknowns each; eliminating them leaves two conditions on $\lambda^{\rm c}$ at each end. At the inlet $g_0=(u(0)-m_*,T(0)-D_4)^T$ , the rows of $g_{0,z}$ are $e_u^T$ and $e_T^T$ , and the last two components of the inlet equation contain no $\eta_0$ :

\begin{gathered}\lambda^{\rm c}_\Pi(0)=-\phi_{0,\Pi},\qquad\lambda^{\rm c}_{\mathcal H}(0)=-\phi_{0,\mathcal H},\\ \eta_0=-\bigl(\lambda^{\rm c}_u(0)+\phi_{0,u},\ \lambda^{\rm c}_T(0)+\phi_{0,T}\bigr)^T.\end{gathered}

加热出口 $g_1=(\Pi-m_*u-D_5,\,q)^T$ 。由 $q=\mathcal H-m_*(c_pT+u^2/2)+u\tau$ 与 $\tau=m_*u+m_*T/u-\Pi$ 逐项求导，再在出口用 $m_*T/u=\rho T=p$ 、 $p-\tau=D_5$ 与 $c_p-1=c_v$ ：

At the heated outlet $g_1=(\Pi-m_*u-D_5,\,q)^T$ . Differentiate $q=\mathcal H-m_*(c_pT+u^2/2)+u\tau$ and $\tau=m_*u+m_*T/u-\Pi$ term by term, then use $m_*T/u=\rho T=p$ , $p-\tau=D_5$ and $c_p-1=c_v$ at the outlet:

\begin{aligned}q_z&=\Bigl(\tau-\frac{m_*T}{u},\ -m_*(c_p-1),\ -u,\ 1\Bigr),\\ q_z(1)&=\bigl(-D_5,\ -m_*c_v,\ -u(1),\ 1\bigr).\end{aligned}

把 $g_{1,z}$ 的两行 $(-m_*,0,1,0)$ 与 $q_z(1)$ 代入 $\lambda^{\rm c}(1)=\phi_{1,z}^T+g_{1,z}^T\eta_1$ ，逐分量写出：

Substitute the two rows of $g_{1,z}$ , $(-m_*,0,1,0)$ and $q_z(1)$ , into $\lambda^{\rm c}(1)=\phi_{1,z}^T+g_{1,z}^T\eta_1$ component by component:

\begin{aligned}\lambda^{\rm c}_{\mathcal H}(1)-\phi_{1,\mathcal H}&=\eta_{1,2},&\lambda^{\rm c}_\Pi(1)-\phi_{1,\Pi}&=\eta_{1,1}-u(1)\,\eta_{1,2},\\ \lambda^{\rm c}_u(1)-\phi_{1,u}&=-m_*\eta_{1,1}-D_5\,\eta_{1,2},&\lambda^{\rm c}_T(1)-\phi_{1,T}&=-m_*c_v\,\eta_{1,2}.\end{aligned}

前两式确定 $\eta_1$ ，代入后两式就是出口的两个伴随条件。所以伴随和原始问题一样，每端两个条件，但给定的是互补的分量：原始问题在入口给定 $u$ 、 $T$ ，伴随在入口给定与自由通量 $\Pi$ 、 $\mathcal H$ 共轭的 $\lambda_\Pi$ 、 $\lambda_{\mathcal H}$ （图 A.1）。制造解算例的出口条件是 $u$ 、 $T$ ，同理出口的伴随条件落在 $\lambda_\Pi$ 、 $\lambda_{\mathcal H}$ 上。

The first two determine $\eta_1$ ; substituted into the last two they are the two outlet adjoint conditions. Like the primal problem, the adjoint thus has two conditions at each end, but on complementary components: the primal prescribes $u$ and $T$ at the inlet, the adjoint prescribes there the components $\lambda_\Pi$ and $\lambda_{\mathcal H}$ conjugate to the free fluxes $\Pi$ and $\mathcal H$ (Figure A.1). In the manufactured case the outlet conditions are on $u$ and $T$ , and likewise the outlet adjoint conditions fall on $\lambda_\Pi$ and $\lambda_{\mathcal H}$ .

图 A.1 · 约化四维问题在两端的条件。实心圆是给定的分量，空心圆是待求的分量。原始问题在入口给定 u、T，自由的 Π(0)、ℋ(0) 由射击法确定，出口满足两个组合条件 g₁ = 0。伴随在入口给定与自由通量共轭的 λ_Π、λ_ℋ，λ_u(0)、λ_T(0) 随后给出 η₀；出口的 λ(1) 含两个待定乘子 η₁，由入口的两个条件解出。两端各两个条件，合计四个，等于 ODE 的维数。

Figure A.1 · The conditions of the reduced four-dimensional problem at each end. Filled circles are prescribed components, open circles unknowns. The primal problem prescribes u and T at the inlet; the free Π(0) and ℋ(0) are found by shooting, and the outlet satisfies two combined conditions g₁ = 0. The adjoint prescribes at the inlet the components λ_Π and λ_ℋ conjugate to the free fluxes, while λ_u(0) and λ_T(0) then give η₀; λ(1) at the outlet contains two undetermined multipliers η₁, solved from the two inlet conditions. Two conditions at each end make four, the dimension of the ODE.

A.4 求端点乘子并保留完整梯度

A.4 Solve endpoint multipliers and retain the full gradient

入口自由变量是两个通量，所以入口条件首先要求 $\lambda^{\rm c}_{\Pi}(0)=-\phi_{0,\Pi}$ 、 $\lambda^{\rm c}_{\mathcal H}(0)=-\phi_{0,\mathcal H}$ 。从出口反向积分一个由目标驱动的特解和两个齐次响应，再用这两条入口条件解 $\eta_1$ 的 $2\times2$ 系统，形成实际终值并重新积分。入口另外两个条件确定 $\eta_0$ 。

The two inlet fluxes are free, so first impose $\lambda^{\rm c}_{\Pi}(0)=-\phi_{0,\Pi}$ and $\lambda^{\rm c}_{\mathcal H}(0)=-\phi_{0,\mathcal H}$ . Integrate one objective-forced particular solution and two homogeneous responses backward from the outlet. Use these two inlet conditions to solve a $2\times2$ system for $\eta_1$ , form the actual terminal value and integrate again. The remaining inlet conditions determine $\eta_0$ .

g^{\rm c}=\phi_{0,D}^T+\phi_{1,D}^T+\int_0^1(\ell_D^T+f_D^T\lambda^{\rm c})\,dx+g_{0,D}^T\eta_0+g_{1,D}^T\eta_1.

入口温度项可化为 $\lambda_T^{\rm c}(0)+\phi_{0,T}$ ；出口参数项必须保留 $g_{1,D}^T\eta_1$ 。入口热流目标会同时改变入口伴随条件与显式目标项。反向积分梯度时积分步长为负，代码最终取负号恢复从 0 到 1 的积分方向，这与任意更改伴随符号无关。

The inlet-temperature contribution reduces to $\lambda_T^{\rm c}(0)+\phi_{0,T}$ ; retain the outlet parameter term $g_{1,D}^T\eta_1$ . An inlet heat-flux objective changes both inlet adjoint conditions and direct objective terms. Gradient quadrature marched backward has negative integration steps, so the implementation negates it to recover the 0-to-1 integral, not to change the adjoint convention arbitrarily.

这是先推导连续 ODE 伴随再作数值积分，不是 RK4 程序的离散转置，也不调用有限体积残差。它仍共用物性及源项定义，独立性体现在方程、边界求解与导数路线。下方可调整 ODE 步数，检查射击残差、连续伴随边界残差及完整重算连续差分，再与正文总设计梯度对照。

This derives a continuous ODE adjoint before numerical integration; it is neither the discrete transpose of RK4 nor a call to the finite-volume residual. Transport and source definitions are shared, while equations, boundary solving and derivative routes are distinct. Vary ODE steps below, inspect shooting and continuous adjoint boundary residuals, and compare fully re-solved continuous differences before comparing total gradients with the main discrete model.

solve the continuous shooting problem for z(x)
compute f_z, f_D and objective-density partials along the trajectory
march one forced adjoint and two homogeneous responses backward
solve the two inlet free-flux conditions for outlet multipliers
march the actual adjoint and gradient quadrature backward
reverse the quadrature sign and add both endpoint contributions
verify by independently re-solving the continuous BVP after perturbation

图 A.2 · 三种目标的连续伴随 λ^c(x)（加热算例，页面默认参数，256 步 RK4）。(a) 温度跟踪由体内源项 −ℓ_z^T 驱动，入口处的 λ_Π、λ_ℋ 为零；(b) 出口温度只由出口条件驱动，体内方程是齐次的；(c) 入口热流在入口给定 λ_Π(0) = m_* = 0.35 与 λ_ℋ(0) = −1。

Figure A.2 · The continuous adjoint λ^c(x) for the three objectives (heating case, page defaults, 256 RK4 steps). (a) Temperature tracking is driven by the interior source −ℓ_z^T, with λ_Π and λ_ℋ zero at the inlet; (b) the outlet temperature is driven only by the outlet condition, the interior equation being homogeneous; (c) the inlet heat flux prescribes λ_Π(0) = m_* = 0.35 and λ_ℋ(0) = −1 at the inlet.

逐步推导：三种目标从哪里驱动伴随

Step by step: where each objective drives the adjoint

三种目标只在 $\ell$ 、 $\phi_0$ 、 $\phi_1$ 中出现一处，伴随方程的其余部分完全相同。温度跟踪 $\ell=(T-T^*)^2/2$ 给出体内源项 $\ell_z=(0,\,T-T^*,\,0,\,0)$ ，两端的 $\phi$ 为零，所以入口 $\lambda^{\rm c}_\Pi(0)=\lambda^{\rm c}_{\mathcal H}(0)=0$ 。出口温度 $\phi_1=T(1)$ 只有 $\phi_{1,z}=e_T^T$ ，体内方程是齐次的，伴随完全由出口条件驱动。入口热流 $\phi_0=q(0)$ 由上面的 $q_z$ 得到：

Each objective appears in exactly one of $\ell$ , $\phi_0$ and $\phi_1$ ; the rest of the adjoint problem is identical. Temperature tracking, $\ell=(T-T^*)^2/2$ , gives the interior source $\ell_z=(0,\,T-T^*,\,0,\,0)$ and zero endpoint terms, so $\lambda^{\rm c}_\Pi(0)=\lambda^{\rm c}_{\mathcal H}(0)=0$ at the inlet. The outlet temperature $\phi_1=T(1)$ has only $\phi_{1,z}=e_T^T$ ; the interior equation is homogeneous and the adjoint is driven entirely by the outlet condition. The inlet heat flux $\phi_0=q(0)$ follows from $q_z$ above:

\lambda^{\rm c}_\Pi(0)=-\frac{\partial q}{\partial\Pi}=u(0)=m_*=0.35,\qquad\lambda^{\rm c}_{\mathcal H}(0)=-\frac{\partial q}{\partial\mathcal H}=-1.

图 A.2(c) 的入口值正是这两个数。入口温度的导数是 A.4 的 $\lambda^{\rm c}_T(0)+\phi_{0,T}$ 。出口温度目标没有 $\phi_0$ ，所以 $dT(1)/dD_4=\lambda^{\rm c}_T(0)=0.9672$ ：入口温度升高 1，出口温度只升高约 0.97。入口热流目标的 $\phi_0$ 还依赖 $u$ 、 $T$ ，它们进入 $\eta_0$ ，不改变上面两个条件。

These are the inlet values in Figure A.2(c). The inlet-temperature derivative is $\lambda^{\rm c}_T(0)+\phi_{0,T}$ from A.4. The outlet-temperature objective has no $\phi_0$ , so $dT(1)/dD_4=\lambda^{\rm c}_T(0)=0.9672$ : raising the inlet temperature by 1 raises the outlet temperature by only about 0.97. The inlet heat-flux objective's $\phi_0$ also depends on $u$ and $T$ ; these enter $\eta_0$ and leave the two conditions above unchanged.

附录 B：黏性通量的局部求导

Appendix B: local viscous-flux differentiation

B.1 对一个实际面检查完整链式法则

B.1 Check the complete chain rule on an actual face

选择面编号后，表中显示该面的离散梯度、物性、应力和热流及其 tangent（代码名 x_d）。当前状态种子固定为 $v_j=\sin(j+1)$ ，净通量的伴随种子固定为 $\overline H_f=(0.3,-0.7,0.4)^T$ ；控件改变面编号，不改变种子。局部反传返回全局状态槽位上的贡献，因为边界和二阶重构可能依赖多个邻点。

Select a face to display its discrete gradients, transport, stress and heat flux with their tangents (code names x_d). The fixed state seed is $v_j=\sin(j+1)$ and the adjoint seed of the net flux is $\overline H_f=(0.3,-0.7,0.4)^T$ . The control changes the face, not the seeds. Local reverse accumulation returns contributions to global state slots because boundaries and second-order reconstruction can depend on several neighbours.

\overline H_f^T\dot H_f=\overline H_f^T(H_{f,U}v)=v^T(H_{f,U}^T\overline H_f).

这个局部点积只检查状态路径，参数种子 $\dot D=0$ 。它不能替代物性参数或热源参数的完整差分检查。边界面也可选，但必须按实际边界闭合解释表中数据；格点端点的梯度是第 4 节的三点单侧差分；随后入口三行与出口动量、能量两行被约束替换，出口质量行仍是平衡行。

This local dot product checks only state paths, with the parameter seed $\dot D=0$ . It does not replace full differences for transport or heat-source parameters. Boundary faces can be selected, but their values follow actual closures; node endpoint gradients use the three-point one-sided differences of Section 4; the three inlet rows and the outlet momentum and energy rows are then replaced by constraints, while the outlet mass row stays a balance row.

图 B.1 · 一个面的通量用到哪些存储状态（格点页）。内部面：一阶只用两侧的 L、R；二阶的对流面状态由重构得到，用到 L−1 到 R+1 四个节点，黏性部分始终只用 L、R。端点面在节点上，梯度用三点单侧差分。入口节点的三行被约束替换，这个面只进入入口热流目标 q(0)；出口面进入出口约束行，使它们依赖 M−3 到 M−1 三个节点。

Figure B.1 · The stored states that one face flux uses (node page). Interior face: first order uses only L and R on either side; second-order convective face states come from reconstruction and use the four nodes L−1 to R+1, while the viscous part always uses L and R only. The end faces sit on nodes and their gradients use three-point one-sided differences. The rows of the inlet node are replaced by constraints, so that face enters only the inlet heat-flux objective q(0); the outlet face enters the outlet constraint rows and makes them depend on the three nodes M−3 to M−1.

逐步推导：面模板决定 $A_h$ 的带宽

Step by step: face stencils fix the bandwidth of $A_h$

第 $i$ 个控制体的残差只含它两侧两个面的通量（再减去源项），所以残差块 $i$ 依赖的状态块，是这两个面模板的并集。内部面的模板取自图 B.1(a)(b)：

The residual of control volume $i$ contains only the fluxes of its two faces (minus sources), so the state blocks that residual block $i$ depends on are the union of the two face stencils. For interior faces, Figure B.1(a)(b) gives:

\begin{aligned}\mathcal S_1(i)&=\{i-1,i\}\cup\{i,i+1\}=\{i-1,i,i+1\},\\ \mathcal S_2(i)&=\{i-2,\ldots,i+1\}\cup\{i-1,\ldots,i+2\}=\{i-2,\ldots,i+2\}.\end{aligned}

按标量计，块 $i$ 的第一行（行号 $3i$ ）最远用到第 $3(i+s)+2$ 列，最后一行（行号 $3i+2$ ）最远用到第 $3(i-s)$ 列，其中一阶 $s=1$ 、二阶 $s=2$ 。所以上、下带宽都是 $3s+2$ ：一阶 5，二阶 8。

Counting scalars, the first row of block $i$ (row $3i$ ) reaches column $3(i+s)+2$ and its last row (row $3i+2$ ) reaches column $3(i-s)$ , with $s=1$ at first order and $s=2$ at second. Both bandwidths are therefore $3s+2$ : 5 at first order and 8 at second.

格点页还要看两端。入口节点的三行被约束 $\rho_0$ 、 $u_0$ 、 $T_0$ 替换，只依赖节点 0。出口加热约束 $p-\tau=D_5$ 、 $q=0$ 用三点单侧差分（图 B.1(d)），依赖节点 $M-3$ 到 $M-1$ ：出口约束行（行号 $3(M-1)+1$ 、 $3(M-1)+2$ ）最远到第 $3(M-3)$ 列，相距 $3\cdot2+2=8$ 。所以一阶矩阵的下带宽是 8 而不是 5，转置后就是附录 C 中 $A_h^T$ 的上带宽 8；二阶时内部已经是 8，不再增加。

The node page also needs the two ends. The inlet node rows are replaced by the constraints on $\rho_0$ , $u_0$ and $T_0$ and depend on node 0 only. The heated outlet constraints $p-\tau=D_5$ and $q=0$ use three-point one-sided differences (Figure B.1(d)) and depend on nodes $M-3$ to $M-1$ : the outlet constraint rows ( $3(M-1)+1$ and $3(M-1)+2$ ) reach column $3(M-3)$ , a distance of $3\cdot2+2=8$ . The first-order matrix therefore has lower bandwidth 8, not 5, which becomes the upper bandwidth 8 of $A_h^T$ in Appendix C; at second order the interior already has 8 and nothing grows.

B.2 热流目标的显式参数导数

B.2 Direct parameter derivatives of the heat-flux objective

J_{q,h}=-\kappa_0G_0^T,\qquad\dot J_{q,h}=-\kappa_0\dot G_0^T-G_0^T\left[\frac{c_p}{Pr}\left(\mu_{T,0}\dot T_0^v+\frac{\mu_0}{\mu_{\rm ref}}\dot\mu_{\rm ref}\right)-\frac{\kappa_0}{Pr}\dot{Pr}\right].

固定状态时，格点梯度与面温度都固定；格心入口面温度和边界梯度仍显式依赖入口温度参数。分别设置状态种子 $v$ 与参数种子 $\dot D$ ，才能辨认 $j_U$ 和 $j_D$ 。总梯度必须同时包含显式目标与残差反传两部分。

At fixed state, node gradients and face temperatures are fixed. At the cell inlet, face temperature and boundary gradient still depend explicitly on the prescribed temperature parameter. Separate the state seed $v$ and the parameter seed $\dot D$ to identify $j_U$ and $j_D$ . The total gradient retains both direct objective and residual-reverse contributions.

B.3 可微范围

B.3 Differentiability scope

本页没有 limiter，开放参数限制在光滑正值分支，物性完整求导。越界或 Newton 失败会清空旧结果并报告失败。若加入截断、裁剪、激波捕捉或分支选择，需要重新定义求导对象和可微范围；当前点积与差分结果不能自动延伸到那些场景。

There is no limiter, parameters target a smooth positive branch, and transport is fully differentiated. Invalid parameters or Newton failure clear stale results and report failure. Clipping, shock capture or branch selection would require a new derivative contract and differentiability analysis; current dot and difference checks do not automatically extend to those settings.

附录 C：显式矩阵与直接伴随求解

Appendix C: explicit matrices and direct adjoint solves

C.1 这里显式组装，不采用 matrix-free 直接求解

C.1 Explicit assembly for the direct solve

(A_h)_{:,j}=A_he_j,\quad j=0,\ldots,3M-1,\qquad \mathsf M=A_h^T,\quad b=-j_U.

每次 JVP 以单位状态种子 $v=e_j$ 、参数种子 $\dot D=0$ 给出 $\dot R_h=A_he_j$ ，即一列。把完整精确矩阵转置后，再作带部分主元的 LU。这里存储了稠密矩阵，所以直接求解路径明确不是 matrix-free；一阶和二阶都使用各自实际残差的精确矩阵，不拿预处理矩阵替代。

Each JVP with a unit state seed $v=e_j$ and the parameter seed $\dot D=0$ gives $\dot R_h=A_he_j$ , one column. Transpose the complete exact matrix, then factor it by partial-pivot LU. This stores a dense matrix and is explicitly not a matrix-free solve. Both spatial orders use their actual residual's exact matrix, not the preconditioning matrix.

\mathsf P_{\rm piv}A_h^T=L_{\rm LU}U_{\rm LU},\qquad L_{\rm LU}y=\mathsf P_{\rm piv}(-j_U),\qquad U_{\rm LU}\psi_h=y.

$\mathsf P_{\rm piv}$ 是行置换矩阵，与 GS 的 $P_h$ 无关； $U_{\rm LU}$ 是上三角因子，与流场状态 $U_i$ 无关。每步选当前列绝对值最大的主元并交换整行，包括之前保存的下三角乘数。随后前代、回代，最后用原始精确转置接口检查真残差。

The row permutation $\mathsf P_{\rm piv}$ is unrelated to the GS matrix $P_h$ , and upper factor $U_{\rm LU}$ is unrelated to flow state $U_i$ . Each stage selects the largest-magnitude pivot in the current column and swaps complete rows, including previously stored lower multipliers. Forward and backward substitution follow, then the original exact transpose interface checks the true residual.

for j in state slots:
    A[:,j] = L.jvp(U_d=unit(j), D_d=0).R_d
M = transpose(A); rhs = -objective_state_gradient
for k = 0..nd-1:
    pivot = argmax(abs(M[k:nd,k]))
    swap complete rows and permutation
    eliminate below pivot; store multipliers
y = forward_substitution(L, permuted_rhs)
psi = backward_substitution(U, y)
check norm(exact_AT(psi) + jU)

小网格适合直接法：组装后可观察每一行，分解后多个目标可复用因子；流场、参数、阶数或边界类型改变就要重新组装。稠密 LU 的存储和运算量分别随未知量数平方、立方增长，不能由本页小模型的速度推断大规模非结构问题也应采用稠密直接法。滑块显示消元阶段；下方解向量表始终显示完成分解后的解，不是滑块阶段的未完成解。

Small meshes suit a direct method: inspect each assembled row and reuse factors for several objectives. Reassemble after changes in flow, parameters, order or boundary type. Dense LU storage and work grow quadratically and cubically with the number of unknowns, so small-model timings do not prescribe dense direct methods for large unstructured problems. The slider shows elimination stages; the solution table always shows the completed solve, not a partial solution at the selected stage.

逐步推导：主元为什么要换行，换行带来多少填充

Step by step: why pivoting swaps rows here, and how much fill it creates

$A_h^T$ 由 $3\times3$ 块组成（本页 N = 12 时 39 维）。附录 B 的面模板给出带宽：一阶下、上带宽为 5、8（上带宽 8 来自出口约束行的三点单侧差分），二阶都是 8。与 Burgers 不同，对角元并不总是本列剩余元素中绝对值最大者，所以部分主元真的会交换行：一阶 12 次、二阶 26 次。一阶第一次换行发生在 $k=8$ ：当前对角元的绝对值是 6.65，第 10 行同列元素的绝对值是 8.15，于是交换这两行。

$A_h^T$ consists of $3\times3$ blocks (dimension 39 here for N = 12). The face stencils of Appendix B give the bandwidths: 5 below and 8 above at first order (the upper bandwidth 8 comes from the three-point one-sided differences of the outlet constraint rows), and 8 on both sides at second order. Unlike Burgers, the diagonal is not always the largest remaining entry of its column, so partial pivoting really swaps rows: 12 times at first order and 26 at second. The first first-order swap happens at $k=8$ : the current diagonal entry has magnitude 6.65 and row 10 has 8.15 in the same column, so the two rows are exchanged.

换行最多让 $U_{\rm LU}$ 变宽多少：带宽为 $p$ 、 $q$ 时，第 $k$ 步的主元只能来自第 $k$ 到 $k+p$ 行（更下面的行在第 $k$ 列为零）。被换到第 $k$ 行的行，非零元至多到第 $k+p+q$ 列，所以 $U_{\rm LU}$ 的上带宽不超过 $p+q$ 。

How much wider can swaps make $U_{\rm LU}$ ? With bandwidths $p$ and $q$ , the pivot of step $k$ can only come from rows $k$ to $k+p$ (rows further down are zero in column $k$ ). A row moved to position $k$ has nonzeros at most up to column $k+p+q$ , so the upper bandwidth of $U_{\rm LU}$ never exceeds $p+q$ .

一阶有行从下面的块行换上来，带着更靠右的非零元： $U_{\rm LU}$ 的上带宽从 8 增加到 9（不超过 $p+q=13$ ），消元还在原来为零的位置产生 29 个新的非零元（填充）。二阶有行从下面的块行换上来，带着更靠右的非零元： $U_{\rm LU}$ 的上带宽从 8 增加到 11（不超过 $p+q=16$ ），消元还在原来为零的位置产生 77 个新的非零元（填充）。图 C.1 是两种阶数的 $A_h^T$ ，图 C.2 是一阶分解后的因子。

At first order rows move up from block rows below, bringing entries further right: the upper bandwidth of $U_{\rm LU}$ grows from 8 to 9 (never beyond $p+q=13$ ), and elimination creates 29 new nonzeros where the matrix had zeros (fill). At second order rows move up from block rows below, bringing entries further right: the upper bandwidth of $U_{\rm LU}$ grows from 8 to 11 (never beyond $p+q=16$ ), and elimination creates 77 new nonzeros where the matrix had zeros (fill). Figure C.1 shows $A_h^T$ for both orders and Figure C.2 the first-order factors.

图 C.1 · 本页格点格式的精确 A_h^T 的非零结构（N = 12，页面默认参数的收敛解），细线划分 3×3 块。

Figure C.1 · The nonzero pattern of the exact A_h^T of this page's node scheme (N = 12, converged state at the page defaults); thin lines mark the 3×3 blocks.

图 C.2 · 一阶 A_h^T 的带主元 LU 因子（按置换后的行顺序）。绿色为 L_LU 的乘数，蓝色为 U_LU；红框是填充，左侧琥珀色圆点标出被交换过的行。

Figure C.2 · The partial-pivot LU factors of the first-order A_h^T (rows in pivoted order). Green are the multipliers of L_LU, blue the entries of U_LU; red outlines mark fill, and amber dots in the margin mark rows moved by swaps.

逐步推导：一次分解同时服务 tangent 与 adjoint

Step by step: one factorisation serves both the tangent and the adjoint

C.1 分解的是 $A_h^T$ 。若改为分解 $A_h$ ，同一组因子也能解伴随：由 $\mathsf P_{\rm piv}A_h=L_{\rm LU}U_{\rm LU}$ 得 $A_h^T=U_{\rm LU}^TL_{\rm LU}^T\mathsf P_{\rm piv}$ ，伴随只是把两次三角求解的顺序和方向换过来：

C.1 factors $A_h^T$ . Had it factored $A_h$ , the same factors would also solve the adjoint: $\mathsf P_{\rm piv}A_h=L_{\rm LU}U_{\rm LU}$ gives $A_h^T=U_{\rm LU}^TL_{\rm LU}^T\mathsf P_{\rm piv}$ , and the adjoint merely swaps the order and direction of the two triangular solves:

\underbrace{L_{\rm LU}y=\mathsf P_{\rm piv}t,\quad U_{\rm LU}v=y}_{A_hv=t=-C_h\dot D}\qquad\qquad\underbrace{U_{\rm LU}^Tz=-j_U,\quad L_{\rm LU}^Ty=z,\quad\psi_h=\mathsf P_{\rm piv}^Ty}_{A_h^T\psi_h=-j_U}.

逐步推导：后向误差换算成前向误差

Step by step: turning the backward error into a forward error

用 C.1 检查的真残差定义后向误差：

Define the backward error from the true residual checked in C.1:

\eta=\frac{\|A_h^T\psi_h+j_U\|_\infty}{\|A_h^T\|_\infty\|\psi_h\|_\infty+\|j_U\|_\infty}.

它恰好是使 $\psi_h$ 成为某个扰动系统精确解的最小相对扰动（Rigal–Gaches）：存在 $E$ 、 $f$ 使 $(A_h^T+E)\psi_h=-j_U+f$ ， $\|E\|_\infty\le\eta\|A_h^T\|_\infty$ ， $\|f\|_\infty\le\eta\|j_U\|_\infty$ 。精确解满足 $A_h^T\psi_h^*=-j_U$ ，两式相减得 $A_h^T(\psi_h-\psi_h^*)=f-E\psi_h$ ，取范数：

It is exactly the smallest relative perturbation for which $\psi_h$ solves a perturbed system (Rigal–Gaches): there are $E$ and $f$ with $(A_h^T+E)\psi_h=-j_U+f$ , $\|E\|_\infty\le\eta\|A_h^T\|_\infty$ and $\|f\|_\infty\le\eta\|j_U\|_\infty$ . The exact solution satisfies $A_h^T\psi_h^*=-j_U$ ; subtracting gives $A_h^T(\psi_h-\psi_h^*)=f-E\psi_h$ , and taking norms:

\frac{\|\psi_h-\psi_h^*\|_\infty}{\|\psi_h\|_\infty}\le\kappa_\infty(A_h^T)\,\eta\Bigl(1+\frac{\|j_U\|_\infty}{\|A_h^T\|_\infty\|\psi_h\|_\infty}\Bigr)\lesssim2\,\kappa_\infty(A_h^T)\,\eta.

本页 N = 12（格点，温度跟踪目标）：一阶 $\kappa_\infty(A_h^T)=2098$ 、 $\eta=2.4\times10^{-17}$ ；二阶 $\kappa_\infty(A_h^T)=2853$ 、 $\eta=4.0\times10^{-17}$ 。伴随的相对前向误差不超过约 $2.3\times10^{-13}$ 。作为检查，用补偿点积算残差、迭代细化三次得到 $\psi_h^*$ ，LU 解的实际相对误差一阶是 $3.2\times10^{-15}$ 、二阶是 $1.1\times10^{-14}$ ，都在界内。

Here N = 12 (node, temperature-tracking objective): first order $\kappa_\infty(A_h^T)=2098$ and $\eta=2.4\times10^{-17}$ ; second order $\kappa_\infty(A_h^T)=2853$ and $\eta=4.0\times10^{-17}$ . The relative forward error of the adjoint is therefore at most about $2.3\times10^{-13}$ . As a check, three steps of iterative refinement with residuals from a compensated dot product give $\psi_h^*$ ; the actual relative error of the LU solution is $3.2\times10^{-15}$ at first order and $1.1\times10^{-14}$ at second, both within the bound.

完整可运行实现Complete runnable implementation

(function(root){'use strict';
const gamma=1.4,alpha=2,zeros=n=>Array(n).fill(0),dot=(a,b)=>a.reduce((s,x,i)=>s+x*b[i],0),norm=a=>Math.sqrt(dot(a,a));
const real={c:x=>x,add:(a,b)=>a+b,sub:(a,b)=>a-b,mul:(a,b)=>a*b,div:(a,b)=>a/b,pow:(a,p)=>a**p,sin:Math.sin,cos:Math.cos};
function tape(){const val=[],pa=[],pb=[],da=[],db=[],inputs=[];function node(v,a=-1,b=-1,x=0,y=0){const i=val.length;val.push(v);pa.push(a);pb.push(b);da.push(x);db.push(y);return i;}const op={c:x=>node(x),input:x=>{const i=node(x);inputs.push(i);return i;},add:(a,b)=>node(val[a]+val[b],a,b,1,1),sub:(a,b)=>node(val[a]-val[b],a,b,1,-1),mul:(a,b)=>node(val[a]*val[b],a,b,val[b],val[a]),div:(a,b)=>node(val[a]/val[b],a,b,1/val[b],-val[a]/val[b]**2),pow:(a,p)=>node(val[a]**p,a,-1,p*val[a]**(p-1)),sin:a=>node(Math.sin(val[a]),a,-1,Math.cos(val[a])),cos:a=>node(Math.cos(val[a]),a,-1,-Math.sin(val[a]))};
function forward(seed){const d=new Float64Array(val.length);inputs.forEach((i,k)=>d[i]=seed[k]||0);for(let i=0;i<val.length;i++){if(pa[i]>=0)d[i]+=da[i]*d[pa[i]];if(pb[i]>=0)d[i]+=db[i]*d[pb[i]];}return d;}
function reverse(outputs,seeds){const w=new Float64Array(val.length);outputs.forEach((i,k)=>w[i]+=seeds[k]);for(let i=val.length-1;i>=0;i--){if(pa[i]>=0)w[pa[i]]+=da[i]*w[i];if(pb[i]>=0)w[pb[i]]+=db[i]*w[i];}return inputs.map(i=>w[i]);}
return {op,val,forward,reverse};}
const complex={c:x=>[x,0],add:(a,b)=>[a[0]+b[0],a[1]+b[1]],sub:(a,b)=>[a[0]-b[0],a[1]-b[1]],mul:(a,b)=>[a[0]*b[0]-a[1]*b[1],a[0]*b[1]+a[1]*b[0]],div:(a,b)=>{const d=b[0]**2+b[1]**2;return [(a[0]*b[0]+a[1]*b[1])/d,(a[1]*b[0]-a[0]*b[1])/d];},pow:(a,p)=>{const r=Math.hypot(...a)**p,t=p*Math.atan2(a[1],a[0]);return [r*Math.cos(t),r*Math.sin(t)];},sin:a=>[Math.sin(a[0])*Math.cosh(a[1]),Math.cos(a[0])*Math.sinh(a[1])],cos:a=>[Math.cos(a[0])*Math.cosh(a[1]),-Math.sin(a[0])*Math.sinh(a[1])]};
function gmres(A,b,P=x=>x.slice(),tol=1e-11,max=b.length){const n=b.length,beta=norm(b);if(beta<=tol)return {x:zeros(n),residual:beta,history:[beta],states:[zeros(n)],converged:true};const V=[b.map(x=>x/beta)],Z=[],H=Array.from({length:max+1},()=>zeros(max)),cs=zeros(max),sn=zeros(max),g=zeros(max+1),history=[beta],estimates=[beta],states=[zeros(n)];g[0]=beta;let m=0;
 for(let j=0;j<max;j++){Z.push(P(V[j]));let w=A(Z[j]);for(let pass=0;pass<2;pass++)for(let i=0;i<=j;i++){const h=dot(V[i],w);H[i][j]+=h;w=w.map((x,k)=>x-h*V[i][k]);}H[j+1][j]=norm(w);const nextnorm=H[j+1][j];if(nextnorm>1e-16)V.push(w.map(x=>x/nextnorm));for(let i=0;i<j;i++){const a=H[i][j],b=H[i+1][j];H[i][j]=cs[i]*a+sn[i]*b;H[i+1][j]=-sn[i]*a+cs[i]*b;}const a=H[j][j],bb=H[j+1][j],r=Math.hypot(a,bb);cs[j]=r?a/r:1;sn[j]=r?bb/r:0;H[j][j]=r;H[j+1][j]=0;g[j+1]=-sn[j]*g[j];g[j]*=cs[j];m=j+1;estimates.push(Math.abs(g[j+1]));const yy=zeros(m);for(let ii=m-1;ii>=0;ii--){let a=g[ii];for(let kk=ii+1;kk<m;kk++)a-=H[ii][kk]*yy[kk];yy[ii]=a/H[ii][ii];}const xx=zeros(n);for(let ii=0;ii<m;ii++)for(let kk=0;kk<n;kk++)xx[kk]+=Z[ii][kk]*yy[ii];states.push(xx);history.push(norm(A(xx).map((v,i)=>v-b[i])));if(history.at(-1)<tol||nextnorm<=1e-16)break;}
 const y=zeros(m);for(let i=m-1;i>=0;i--){let a=g[i];for(let j=i+1;j<m;j++)a-=H[i][j]*y[j];y[i]=a/H[i][i];}const x=zeros(n);for(let i=0;i<m;i++)for(let j=0;j<n;j++)x[j]+=Z[i][j]*y[i];const residual=norm(A(x).map((z,i)=>z-b[i]));return {x,residual,history,estimates,states,converged:residual<=Math.max(tol*5,1e-13)};}
const api={real,complex,tape,gmres,dot,norm};if(typeof module!=='undefined')module.exports=api;else root.QMath=api;
})(typeof window!=='undefined'?window:globalThis);

/* Small educational systems: matrix-free exact products; explicit matrices only for LU/conditioning diagnostics. */
(function(root){'use strict';
const Q=typeof module!=='undefined'?require('../q1d/reconstruction/math'):root.QMath;
const {norm,dot}=Q,zeros=n=>Array(n).fill(0),transpose=A=>A[0].map((_,j)=>A.map(r=>r[j])),mv=(A,x)=>A.map(r=>dot(r,x));
function factor(A){const n=A.length,lu=A.map(r=>r.slice()),p=Array.from({length:n},(_,i)=>i),scale=Math.max(...A.flat().map(Math.abs));for(let k=0;k<n;k++){let j=k;for(let i=k+1;i<n;i++)if(Math.abs(lu[i][k])>Math.abs(lu[j][k]))j=i;if(!Number.isFinite(lu[j][k])||Math.abs(lu[j][k])<1e-14*scale)throw Error('Singular LU pivot '+k);[lu[j],lu[k]]=[lu[k],lu[j]];[p[j],p[k]]=[p[k],p[j]];for(let i=k+1;i<n;i++){lu[i][k]/=lu[k][k];for(let j=k+1;j<n;j++)lu[i][j]-=lu[i][k]*lu[k][j];}}return {lu,p};}
function solve(F,b){const {lu,p}=F,n=b.length,x=p.map(i=>b[i]);for(let i=0;i<n;i++)for(let j=0;j<i;j++)x[i]-=lu[i][j]*x[j];for(let i=n-1;i>=0;i--){for(let j=i+1;j<n;j++)x[i]-=lu[i][j]*x[j];x[i]/=lu[i][i];}return x;}
function assemble(L,n){const A=Array.from({length:n},()=>zeros(n));for(let j=0;j<n;j++){const e=zeros(n);e[j]=1;L.jvp(e).R.forEach((v,i)=>A[i][j]=v);}return A;}
function limiter(a,b,kind='minmod',eps=.01){if(kind==='first')return {s:0,da:0,db:0,branch:'first',margin:1};if(kind==='none')return {s:(a+b)/2,da:.5,db:.5,branch:'unlimited',margin:1};if(kind==='smooth'){const den=a*a+b*b+eps*eps,num=a*b*(a+b);return {s:num/den,da:((2*a*b+b*b)*den-2*a*num)/(den*den),db:((a*a+2*a*b)*den-2*b*num)/(den*den),branch:'smooth',margin:1};}const margin=Math.min(Math.abs(a),Math.abs(b),Math.abs(Math.abs(a)-Math.abs(b)));if(a*b<=0)return {s:0,da:0,db:0,branch:'zero',margin};return Math.abs(a)<=Math.abs(b)?{s:a,da:1,db:0,branch:'left',margin}:{s:b,da:0,db:1,branch:'right',margin};}
function slope(o,a,b,kind,eps,mode='full'){const r=limiter(o.value(a),o.value(b),kind,eps),{c,add,mul,sub,div}=o;if(mode==='frozen'){const d=(o.value(a)+o.value(b))/2,phi=Math.abs(d)>1e-15?r.s/d:0;return add(c(r.s),mul(c(phi/2),add(sub(a,c(o.value(a))),sub(b,c(o.value(b))))));}if(kind==='smooth'){const den=add(add(mul(a,a),mul(b,b)),c(eps*eps));return div(mul(mul(a,b),add(a,b)),den);}return add(c(r.s),add(mul(c(r.da),sub(a,c(o.value(a)))),mul(c(r.db),sub(b,c(o.value(b))))));}
function operators(U,p,evalFn,mode='full'){const t=Q.tape();t.op.value=i=>t.val[i];const ui=U.map(t.op.input),pi=p.map(t.op.input),r=evalFn(t.op,ui,pi,mode);return {R:r.R.map(i=>t.val[i]),J:t.val[r.J],jvp:(v,dp=zeros(p.length))=>{const z=t.forward([...v,...dp]);return {R:r.R.map(i=>z[i]),J:z[r.J]};},vjp:(w,j=0)=>{const z=t.reverse([...r.R,r.J],[...w,j]);return {U:z.slice(0,U.length),p:z.slice(U.length)};}};}
const real={...Q.real,value:x=>x};
function gsMatrix(P,sweeps=2){const n=P.length,B=Array.from({length:n},()=>zeros(n));for(let k=0;k<n;k++){const b=zeros(n),x=zeros(n);b[k]=1;for(let s=0;s<sweeps;s++)for(let i=0;i<n;i++){let v=b[i];for(let j=0;j<n;j++)if(i!==j)v-=P[i][j]*x[j];if(Math.abs(P[i][i])<1e-14)throw Error('Zero GS diagonal');x[i]=v/P[i][i];}x.forEach((v,i)=>B[i][k]=v);}return B;}
function compare(L,P,rhs,adj=false,omega=.7,preAction=null){const n=rhs.length,A=assemble(L,n),AA=adj?transpose(A):A,F=factor(AA),exact=solve(F,rhs),B0=preAction?transpose(Array.from({length:n},(_,j)=>{const e=zeros(n);e[j]=1;return preAction(e);})):gsMatrix(P),B=adj?transpose(B0):B0,op=v=>adj?L.vjp(v).U:L.jvp(v).R,pre=v=>mv(B,v),gm=Q.gmres(op,rhs,pre,1e-11,n);let x=zeros(n),history=[],reason='limit';for(let k=0;k<250;k++){const ax=op(x),r=rhs.map((b,i)=>b-ax[i]),rn=norm(r);history.push(rn);if(rn<1e-10){reason='converged';break;}if(!Number.isFinite(rn)||rn>Math.max(1,norm(rhs))*1e8){reason='growth';break;}const d=pre(r);x=x.map((v,i)=>v+omega*d[i]);}
const inv=transpose(Array.from({length:n},(_,j)=>{const e=zeros(n);e[j]=1;return solve(F,e);})),ni=M=>Math.max(...M.map(r=>r.reduce((s,v)=>s+Math.abs(v),0))),cond=ni(AA)*ni(inv),M=AA.map((_,i)=>AA[0].map((_,j)=>(i===j?1:0)-omega*B[i].reduce((s,v,k)=>s+v*AA[k][j],0)));return {exact,luResidual:norm(op(exact).map((v,i)=>v-rhs[i])),gm,gs:{x,history,reason},conditionInfinity:cond,iterationNormInfinity:ni(M),A};}
const api={...Q,zeros,transpose,mv,factor,solve,assemble,limiter,slope,operators,real,gsMatrix,compare};if(typeof module!=='undefined')module.exports=api;else root.EduNumerics=api;
})(typeof window!=='undefined'?window:globalThis);

/* Plane 1D laminar compressible Navier–Stokes. R = Fc_right-Fc_left-Fv_right+Fv_left-S. */
(function(root){'use strict';
const Q=typeof module!=='undefined'?require('../q1d/reconstruction/math'):root.QMath;
const N=typeof module!=='undefined'?require('../experiments/numerics'):root.EduNumerics;
const {norm,dot,real,complex,tape,gmres}=Q,{zeros,assemble,transpose,factor,solve:lu,mv}=N;
const cp=3.5,cv=2.5,mass=.35;
function physics(o,variable){const {c,add,sub,mul,div,pow}=o;
 const prim=U=>{const rho=U[0],u=div(U[1],rho),T=div(sub(div(U[2],rho),mul(c(.5),mul(u,u))),c(cv));return {rho,u,T,p:mul(rho,T)};};
 const cons=(rho,u,T)=>[rho,mul(rho,u),mul(rho,add(mul(c(cv),T),mul(c(.5),mul(u,u))))];
 const transport=(T,b)=>{const mu=variable?mul(b[2],div(mul(c(1.4),pow(T,1.5)),add(T,c(.4)))):b[2];return {mu,kappa:div(mul(c(cp),mu),b[3])};};
 const flux=U=>{const q=prim(U);return [U[1],add(mul(U[1],q.u),q.p),mul(q.u,add(U[2],q.p))];};
 return {prim,cons,transport,flux};
}
function exact(x){const u=mass+.03*Math.sin(Math.PI*x)+.01*Math.sin(2*Math.PI*x),T=1+.05*Math.sin(Math.PI*x)+.02*Math.sin(2*Math.PI*x);return [mass/u,mass,(mass/u)*(cv*T+.5*u*u)];}
function manufactured(x,mu0,Pr,variable){const s=Math.sin(Math.PI*x),c=Math.cos(Math.PI*x),u=mass+.03*s+.01*Math.sin(2*Math.PI*x),T=1+.05*s+.02*Math.sin(2*Math.PI*x),ux=.03*Math.PI*c+.02*Math.PI*Math.cos(2*Math.PI*x),tx=.05*Math.PI*c+.04*Math.PI*Math.cos(2*Math.PI*x),uxx=-.03*Math.PI**2*s-.04*Math.PI**2*Math.sin(2*Math.PI*x),txx=-.05*Math.PI**2*s-.08*Math.PI**2*Math.sin(2*Math.PI*x);
 const mu=variable?mu0*1.4*T**1.5/(T+.4):mu0,mut=variable?mu*(1.5/T-1/(T+.4)):0,k=cp*mu/Pr,kt=cp*mut/Pr,tau=4/3*mu*ux,taux=4/3*(mut*tx*ux+mu*uxx);
 return [0,mass*ux+mass*(tx/u-T*ux/u**2)-taux,mass*(cp*tx+u*ux)-ux*tau-u*taux-kt*tx*tx-k*txx];}
function model({node=false,n=12,order=2,variable=false,mode='heat',objective='tracking',mu=.15,Pr=.72}={}){
 if(n<4)throw Error('At least four intervals required');const h=1/n,m=node?n+1:n,nd=3*m,x=Array.from({length:m},(_,i)=>(i+(node?0:.5))*h),fx=node?[0,...Array.from({length:n},(_,i)=>(i+.5)*h),1]:Array.from({length:n+1},(_,i)=>i*h),vol=x.map((_,i)=>fx[i+1]-fx[i]);
 const params=[mode==='heat'?.03:0,0,mu,Pr,1,1],names=['a₀','a₁','μ_ref','Pr','T_in',mode==='heat'?'p_b':'T_out'];
 function source(o,xx,b){const {c,add,mul}=o,base=mode==='mms'?manufactured(xx,mu,Pr,variable):[0,0,0];return [c(0),c(base[1]),add(c(base[2]),add(mul(b[0],c(Math.sin(Math.PI*xx)**2)),mul(b[1],c(Math.sin(2*Math.PI*xx)**2))))];}
 function weights(i,xx,q=order){if(q===1)return [[i,1]];const d=(xx-x[i])/h;if(i===0)return [[0,1-d],[1,d]];if(i===m-1)return [[i,1+d],[i-1,-d]];return [[i,1],[i-1,-d/2],[i+1,d/2]];}
 function evaluate(o,U,b=params,q=order){const {c,add,sub,mul,div,pow}=o,ph=physics(o,variable),R=Array.from({length:nd},()=>c(0)),local=[],sources=[],qstates=Array.from({length:m},(_,i)=>ph.prim(U.slice(3*i,3*i+3)));
 const sum=a=>a.reduce((s,v)=>add(s,v),c(0)),lin=(a,w)=>sum(w.map(([i,v])=>mul(c(v),a[i]))),rec=(i,xx)=>[0,1,2].map(k=>lin(U.filter((_,j)=>j%3===k),weights(i,xx,q))),at=(key,i)=>qstates[i][key];
 const inlet=ph.cons(c(1),c(mass),b[4]);
 for(let f=0;f<=m;f++){const L=f-1,Ri=f,xx=fx[f];let ul,ur,face,ux,tx;
  if(f>0&&f<m){ul=rec(L,xx);ur=rec(Ri,xx);face=ph.prim([0,1,2].map(k=>mul(c(.5),add(U[3*L+k],U[3*Ri+k]))));ux=mul(c(1/h),sub(at('u',Ri),at('u',L)));tx=mul(c(1/h),sub(at('T',Ri),at('T',L)));}
  else if(node){const i=f===0?0:m-1;ul=ur=U.slice(3*i,3*i+3);face=ph.prim(ul);const w=f===0?[[0,-1.5/h],[1,2/h],[2,-.5/h]]:[[m-1,1.5/h],[m-2,-2/h],[m-3,.5/h]];ux=lin(qstates.map(s=>s.u),w);tx=lin(qstates.map(s=>s.T),w);}
  else if(f===0){ul=ur=inlet;face=ph.prim(inlet);ux=mul(c(1/h),sum([mul(c(-8/3),face.u),mul(c(3),at('u',0)),mul(c(-1/3),at('u',1))]));tx=mul(c(1/h),sum([mul(c(-8/3),face.T),mul(c(3),at('T',0)),mul(c(-1/3),at('T',1))]));}
  else {const ext=key=>sub(mul(c(1.5),at(key,m-1)),mul(c(.5),at(key,m-2)));const rho=ext('rho'),u=mode==='mms'?c(mass):ext('u'),T=mode==='mms'?b[5]:ext('T');ul=ur=ph.cons(rho,u,T);face=ph.prim(ul);ux=mode==='mms'?mul(c(1/h),sum([mul(c(8/3),u),mul(c(-3),at('u',m-1)),mul(c(1/3),at('u',m-2))])):c(0);tx=mode==='mms'?mul(c(1/h),sum([mul(c(8/3),T),mul(c(-3),at('T',m-1)),mul(c(1/3),at('T',m-2))])):c(0);}
  const tr=ph.transport(face.T,b);let tau=mul(c(4/3),mul(tr.mu,ux)),heat=mul(c(-1),mul(tr.kappa,tx));
  if(!node&&f===m&&mode==='heat'){tau=sub(face.p,b[5]);ux=div(tau,mul(c(4/3),tr.mu));}
  let fc=ph.flux(ul),alpha=c(0);if(f>0&&f<m){const a=ph.prim(ul),bb=ph.prim(ur),fr=ph.flux(ur);alpha=mul(c(.5),sum([a.u,bb.u,pow(mul(c(1.4),a.T),.5),pow(mul(c(1.4),bb.T),.5)]));fc=fc.map((v,k)=>sub(mul(c(.5),add(v,fr[k])),mul(mul(c(.5),alpha),sub(ur[k],ul[k]))));}
  const fv=[c(0),tau,sub(mul(face.u,tau),heat)],H=fc.map((v,k)=>sub(v,fv[k]));
  H.forEach((v,k)=>{if(L>=0)R[3*L+k]=add(R[3*L+k],v);if(Ri<m)R[3*Ri+k]=sub(R[3*Ri+k],v);});local.push({x:xx,ul,ur,fc,fv,H,ux,tx,tau,heat,mu:tr.mu,kappa:tr.kappa,u:face.u,T:face.T,p:face.p,alpha});
 }
 let J=c(0);for(let i=0;i<m;i++){const mid=(fx[i]+fx[i+1])/2,half=vol[i]/2,S=[c(0),c(0),c(0)];for(const z of [-1/Math.sqrt(3),1/Math.sqrt(3)]){const xx=mid+half*z,ss=source(o,xx,b);ss.forEach((v,k)=>S[k]=add(S[k],mul(c(half),v)));const T=ph.prim(rec(i,xx)).T,d=sub(T,c(1+.04*Math.sin(Math.PI*xx)));if(objective==='tracking')J=add(J,mul(c(.5*half),mul(d,d)));}S.forEach((v,k)=>R[3*i+k]=sub(R[3*i+k],v));sources.push(S);}
 if(objective==='outlet')J=local[m].T;if(objective==='heatflux')J=local[0].heat;
 const raw=R.slice(),constraints=[];if(node){const put=(row,v)=>{R[row]=v;constraints.push(row);};put(0,sub(qstates[0].rho,c(1)));put(1,sub(qstates[0].u,c(mass)));put(2,sub(qstates[0].T,b[4]));if(mode==='mms'){put(3*(m-1)+1,sub(qstates[m-1].u,c(mass)));put(3*(m-1)+2,sub(qstates[m-1].T,b[5]));}else{put(3*(m-1)+1,sub(sub(local[m].p,local[m].tau),b[5]));put(3*(m-1)+2,local[m].heat);}}
 return {R,J,raw,local,sources,constraints};
 }
 function linearize(U,b=params,q=order){const t=tape(),u=U.map(t.op.input),p=b.map(t.op.input),out=evaluate(t.op,u,p,q),vals=a=>a.map(i=>t.val[i]);
 const extract=(a,d)=>Object.fromEntries(Object.entries(a).map(([k,v])=>[k,k==='x'?v:Array.isArray(v)?v.map(i=>d[i]):d[v]]));
 /* Tapenade-style names: x_d is the tangent (dot) of x, x_b its adjoint (bar). jvp(U_d,D_d) returns the tangents
    R_d, J_d, ... under the keys R, raw, J, local, sources; vjp(R_b,J_b) returns U_b and D_b under the keys U and p. */
 return {R:vals(out.R),J:t.val[out.J],jvp:(U_d,D_d=zeros(b.length))=>{const x_d=t.forward([...U_d,...D_d]);return {R:out.R.map(i=>x_d[i]),raw:out.raw.map(i=>x_d[i]),J:x_d[out.J],local:out.local.map(a=>extract(a,x_d)),sources:out.sources.map(a=>a.map(i=>x_d[i]))};},vjp:(R_b,J_b=0)=>{const x_b=t.reverse([...out.R,out.J],[...R_b,J_b]);return {U:x_b.slice(0,nd),p:x_b.slice(nd)};},constraintVjp:R_b=>t.reverse(out.constraints.map(i=>out.R[i]),out.constraints.map(i=>R_b[i])).slice(0,nd),faceVjp:(f,H_b)=>t.reverse(out.local[f].H,H_b).slice(0,nd)};
 }
 function precondition(U,b=params){const P=assemble(linearize(U,b,1),nd),blocks=Array.from({length:m},(_,i)=>factor(P.slice(3*i,3*i+3).map(r=>r.slice(3*i,3*i+3))));
 // One forward block GS sweep B=(diag(P)+L)^-1; reverse sweep is exactly B^T.
 const apply=(r,tr=false)=>{const z=zeros(nd);for(let it=0;it<m;it++){const i=tr?m-1-it:it,rr=r.slice(3*i,3*i+3);for(let j=tr?i+1:0;tr?j<m:j<i;j++)for(let k=0;k<3;k++)for(let l=0;l<3;l++)rr[k]-=(tr?P[3*j+l][3*i+k]:P[3*i+k][3*j+l])*z[3*j+l];const block=tr?factor(transpose(P.slice(3*i,3*i+3).map(r=>r.slice(3*i,3*i+3)))):blocks[i],v=lu(block,rr);v.forEach((a,k)=>z[3*i+k]=a);}return z;};return {P,apply};}
 function linearSolve(L,rhs,pre,trans=false,method='gmres',tol=1e-11){const op=v=>trans?L.vjp(v).U:L.jvp(v).R;
 if(method==='lu'){const A=assemble(L,nd),x=lu(factor(trans?transpose(A):A),rhs);return {x,history:[norm(rhs),norm(op(x).map((v,i)=>v-rhs[i]))],converged:true,A};}
 if(method==='gmres')return gmres(op,rhs,r=>pre.apply(r,trans),tol,nd);
 let x=zeros(nd),history=[];for(let it=0;it<1500;it++){const ax=op(x),r=rhs.map((v,i)=>v-ax[i]),rn=norm(r);history.push(rn);if(rn<tol)return {x,history,converged:true};if(!Number.isFinite(rn)||rn>1e8)break;const d=pre.apply(r,trans);x=x.map((v,i)=>v+.65*d[i]);}return {x,history,converged:false};}
 function valid(U,b=params){try{const z=evaluate(real,U,b);return z.local.every(f=>[f.ul,f.ur].every(v=>{const s=physics(real,variable).prim(v);return s.rho>.1&&s.T>.3&&s.u>.05&&s.u<.9*Math.sqrt(1.4*s.T);}));}catch{return false;}}
 function initial(){return x.flatMap(xx=>mode==='mms'?exact(xx):[1,mass,cv+.5*mass**2]);}
 function solve(b=params,start,method='gmres',tol=2e-12){let U=start?start.slice():initial();const history=[],states=[];for(let it=0;it<35;it++){const L=linearize(U,b),rn=norm(L.R);history.push(rn);states.push(U.slice());if(rn<tol)return {U,history,states,converged:true};const pre=precondition(U,b),v=linearSolve(L,L.R.map(x=>-x),pre,false,method,2e-13);if(!v.converged)throw Error('Primal linear solve failed');let step=1,ok=false;while(step>1e-8){const trial=U.map((x,i)=>x+step*v.x[i]);if(valid(trial,b)&&norm(evaluate(real,trial,b).R)<(1-1e-4*step)*rn){U=trial;ok=true;break;}step*=.5;}if(!ok)throw Error('No admissible Newton correction');}throw Error('Primal residual did not converge');}
 function sensitivity(sol,b=params,k=0,method='gmres'){const L=linearize(sol.U,b),pre=precondition(sol.U,b),D_d=zeros(b.length);D_d[k]=1;const tangent=linearSolve(L,L.jvp(zeros(nd),D_d).R.map(x=>-x),pre,false,method),adjoint=linearSolve(L,L.vjp(zeros(nd),1).U.map(x=>-x),pre,true,method);if(!tangent.converged||!adjoint.converged)throw Error('Derivative solve did not converge');return {L,pre,tangent,adjoint,gradient:L.vjp(adjoint.x,1).p,directional:L.jvp(tangent.x,D_d).J};}
 function complexGradient(sol,b=params,k=0){const h=1e-25,p=b.map((v,i)=>[v,i===k?h:0]),U=sol.U.map(v=>[v,0]),L=linearize(sol.U,b),pre=precondition(sol.U,b);for(let it=0;it<4;it++){const z=evaluate(complex,U,p),r=z.R.map(v=>v[1]/h);if(norm(r)<1e-10)return {gradient:z.J[1]/h,residual:norm(r)};const s=linearSolve(L,r.map(v=>-v),pre,false,'gmres',1e-13);if(!s.converged)throw Error('Complex correction failed');U.forEach((v,i)=>v[1]+=h*s.x[i]);}const z=evaluate(complex,U,p),residual=norm(z.R.map(v=>v[1]/h));if(residual>1e-9)throw Error('Complex residual not converged');return {gradient:z.J[1]/h,residual};}
 return {node,n,m,nd,h,x,fx,vol,order,mode,variable,objective,params,names,source,weights,initial,valid,solve,linearize,precondition,linearSolve,sensitivity,complexGradient,evaluate:(U,b=params)=>evaluate(real,U,b),evaluateWith:evaluate,physics:physics(real,variable)};
}
const api={model,physics,exact,manufactured,cp,cv,mass};if(typeof module!=='undefined')module.exports=api;else root.NS1D=api;
})(typeof window!=='undefined'?window:globalThis);

/* Independent spatial ODE BVP and continuous ODE adjoint. No FV residual/Jacobian calls. */
(function(root){'use strict';const Q=typeof module!=='undefined'?require('../q1d/reconstruction/math'):root.QMath,N=typeof module!=='undefined'?require('../experiments/numerics'):root.EduNumerics,C=typeof module!=='undefined'?require('./core'):root.NS1D;
const {real,tape,norm,dot}=Q,{zeros,factor,solve:lu}=N,m=C.mass,cp=C.cp;
const add=(a,b,s=1)=>a.map((v,i)=>v+s*b[i]);
function rk(f,x,y,h){const a=f(x,y),b=f(x+h/2,add(y,a,h/2)),c=f(x+h/2,add(y,b,h/2)),d=f(x+h,add(y,c,h));return y.map((v,i)=>v+h*(a[i]+2*b[i]+2*c[i]+d[i])/6);}
function equations(model,o,x,z,b){const {c,add,sub,mul,div}=o,[u,T,P,H]=z,tr=C.physics(o,model.variable).transport(T,b),tau=sub(add(mul(c(m),u),div(mul(c(m),T),u)),P),ux=div(tau,mul(c(4/3),tr.mu)),ht=mul(c(m),add(mul(c(cp),T),mul(c(.5),mul(u,u)))),heat=add(sub(H,ht),mul(u,tau)),tx=div(mul(c(-1),heat),tr.kappa),S=model.source(o,x,b),d=sub(T,c(1+.04*Math.sin(Math.PI*x))),j=model.objective==='tracking'?mul(c(.5),mul(d,d)):c(0);return {f:[ux,tx,S[1],S[2]],j,heat,tau};}
function boundary(model,o,z,b){const {c,sub,mul}=o;return model.mode==='mms'?[sub(z[0],c(m)),sub(z[1],b[5])]:[sub(sub(z[2],mul(c(m),z[0])),b[5]),equations(model,o,1,z,b).heat];}
function derivatives(model,x,z,b){const t=tape(),zz=z.map(t.op.input),bb=b.map(t.op.input),out=equations(model,t.op,x,zz,bb),ids=[...out.f,out.j],jac=ids.map(i=>t.reverse([i],[1]));return {f:out.f.map(i=>t.val[i]),j:t.val[out.j],A:jac.slice(0,4).map(r=>r.slice(0,4)),B:jac.slice(0,4).map(r=>r.slice(4)),jz:jac[4].slice(0,4),jp:jac[4].slice(4)};}
function endpoint(model,z,b,which){const t=tape(),zz=z.map(t.op.input),bb=b.map(t.op.input),phi=which===0&&model.objective==='heatflux'?equations(model,t.op,0,zz,bb).heat:which===1&&model.objective==='outlet'?zz[1]:t.op.c(0),g=boundary(model,t.op,zz,bb),d=t.reverse([phi],[1]);return {phi:t.val[phi],phiz:d.slice(0,4),phip:d.slice(4),g:g.map(i=>t.val[i]),G:g.map(i=>t.reverse([i],[1]).slice(0,4)),Gp:g.map(i=>t.reverse([i],[1]).slice(4))};}
function integrate(model,b,initial,steps,variations=false){let y=[m,b[4],...initial];if(variations)y.push(...[0,0,0,0,1,0,0,1]);const states=[y.slice(0,4)],f=(x,y)=>{const d=derivatives(model,x,y.slice(0,4),b);return variations?[...d.f,...d.A.flatMap(r=>[0,1].map(k=>r.reduce((s,v,j)=>s+v*y[4+2*j+k],0)))]:d.f;};for(let i=0;i<steps;i++){y=rk(f,i/steps,y,1/steps);if(!y.every(Number.isFinite)||y[0]<=.02||y[1]<=.2||y[0]>.95*Math.sqrt(1.4*y[1]))throw Error('Continuous trajectory outside positive subsonic branch');states.push(y.slice(0,4));}return {states,Y:variations?Array.from({length:4},(_,i)=>y.slice(4+2*i,6+2*i)):null};}
function solve(model,b=model.params,steps=256,warm=null){if(!warm&&model.mode==='heat'&&Math.max(Math.abs(b[0]),Math.abs(b[1]))>.04){const base=b.slice();base[0]=base[1]=0;let r=solve(model,base,steps);for(let k=1;k<=6;k++){const next=b.slice();next[0]*=k/6;next[1]*=k/6;r=solve(model,next,steps,r.states[0].slice(2));}return r;}let guess=[m*m+1,m*(cp*b[4]+.5*m*m)];if(model.mode==='mms'){const z=C.exact(0),q=C.physics(real,model.variable).prim(z),mu=C.physics(real,model.variable).transport(q.T,b);const tau=4/3*mu.mu*.05*Math.PI;guess=[m*m+1-tau,m*(cp+.5*m*m)-m*tau-mu.kappa*.09*Math.PI];}
 if(warm)guess=warm.slice();let run,history=[];for(let it=0;it<20;it++){run=integrate(model,b,guess,steps,true);const ep=endpoint(model,run.states.at(-1),b,1),rn=norm(ep.g);history.push(rn);if(rn<2e-12)break;const A=ep.G.map(r=>[0,1].map(k=>r.reduce((s,v,j)=>s+v*run.Y[j][k],0))),d=lu(factor(A),ep.g.map(v=>-v));let step=1,ok=false;while(step>1e-7){try{const trial=add(guess,d,step),r=integrate(model,b,trial,steps),err=norm(endpoint(model,r.states.at(-1),b,1).g);if(err<rn){guess=trial;ok=true;break;}}catch{}step*=.5;}if(!ok)throw Error('Continuous shooting failed');}
 if(history.at(-1)>2e-11)throw Error('Continuous boundary residual not converged');
 const states=run.states,fs=states.map((z,i)=>equations(model,real,i/steps,z,b).f);
 function at(x){const i=Math.min(steps-1,Math.max(0,Math.floor(x*steps))),t=x*steps-i,h=1/steps;return states[i].map((v,k)=>(2*t**3-3*t*t+1)*v+(t**3-2*t*t+t)*h*fs[i][k]+(-2*t**3+3*t*t)*states[i+1][k]+(t**3-t*t)*h*fs[i+1][k]);}
 let integral=0;for(let i=0;i<steps;i++){const x=i/steps,h=1/steps;integral+=h/6*(equations(model,real,x,states[i],b).j+4*equations(model,real,x+h/2,at(x+h/2),b).j+equations(model,real,x+h,states[i+1],b).j);}
 return {states,at,steps,history,J:integral+endpoint(model,states[0],b,0).phi+endpoint(model,states.at(-1),b,1).phi};}
function adjoint(model,ref,b=model.params){const n=ref.steps,end=endpoint(model,ref.states.at(-1),b,1),start=endpoint(model,ref.states[0],b,0);
 const march=(terminal,forcing=true,withGradient=false)=>{let y=[...terminal,...(withGradient?zeros(b.length):[])],values=Array(n+1);values[n]=terminal;const f=(x,y)=>{const d=derivatives(model,x,ref.at(x),b),v=Array.from({length:4},(_,i)=>-d.A.reduce((s,r,j)=>s+r[i]*y[j],0)-(forcing?d.jz[i]:0));if(withGradient)v.push(...b.map((_,i)=>d.jp[i]+d.B.reduce((s,r,j)=>s+r[i]*y[j],0)));return v;};for(let i=n;i>0;i--){y=rk(f,i/n,y,-1/n);values[i-1]=y.slice(0,4);}return {y,values};};
 const a=march(end.phiz),r=end.G.map(v=>march(v,false)),rhs=[2,3].map(i=>-start.phiz[i]-a.y[i]),mat=[2,3].map(i=>r.map(v=>v.y[i])),eta=lu(factor(mat),rhs),terminal=end.phiz.map((v,i)=>v+end.G.reduce((s,r,j)=>s+eta[j]*r[i],0)),out=march(terminal,true,true),l0=out.y.slice(0,4),gradient=b.map((_,i)=>start.phip[i]+end.phip[i]-out.y[4+i]+eta.reduce((s,v,j)=>s+v*end.Gp[j][i],0)+(i===4?l0[1]+start.phiz[1]:0));
 return {values:out.values,eta,gradient,boundaryResidual:Math.hypot(l0[2]+start.phiz[2],l0[3]+start.phiz[3])};}
const api={solve,adjoint,equations,derivatives,boundary,endpoint,rk};if(typeof module!=='undefined')module.exports=api;else root.NSContinuous=api;
})(typeof window!=='undefined'?window:globalThis);

(function(root){'use strict';const C=typeof module!=='undefined'?require('./core'):root.NS1D,N=typeof module!=='undefined'?require('../experiments/numerics'):root.EduNumerics;
function reference(m,i){if(m.node)return C.exact(m.x[i]);const a=m.fx[i],b=m.fx[i+1],r=[0,0,0],z=[-.9602898564975363,-.7966664774136267,-.525532409916329,-.1834346424956498,.1834346424956498,.525532409916329,.7966664774136267,.9602898564975363],w=[.1012285362903763,.2223810344533745,.3137066458778873,.362683783378362,.362683783378362,.3137066458778873,.2223810344533745,.1012285362903763];z.forEach((v,j)=>C.exact((a+b)/2+(b-a)*v/2).forEach((s,k)=>r[k]+=w[j]*s/2));return r;}
function accuracy(options,grids=[8,16,32,64]){const rows=[];for(const order of [1,2]){let prev;for(const n of grids){const m=C.model({...options,mode:'mms',n,order}),s=m.solve(),e=Math.sqrt(m.x.reduce((ss,_,i)=>ss+m.vol[i]*reference(m,i).reduce((sum,v,k)=>sum+(v-s.U[3*i+k])**2,0),0));rows.push({order,n,error:e,rate:prev?Math.log2(prev/e):null,residual:s.history.at(-1)});prev=e;}}return rows;}
/* The players' steps. Tapenade-style labels: x_d is the tangent (dot) of x, x_b its adjoint (bar); the tangent
   player seeds U_d = v, D_d = 0, the adjoint player R_b with J_b = 0. */
function trace(m,sol,kind){const U=sol.U,L=m.linearize(U),raw=m.evaluate(U),U_d=U.map((_,i)=>Math.sin(i+1)),R_b=U.map((_,i)=>Math.cos(2*i+1)),tan=L.jvp(U_d),out=N.zeros(m.nd),steps=[];const step=(zh,en,face,details,updates=[])=>{for(const [i,value,replace]of updates)out[i]=replace?value:out[i]+value;steps.push({zh,en,face,details,values:out.slice()});},nm=s=>kind==='tangent'?s.replace(/^[A-Za-z]+/,a=>a+'_d'):s;
 step('清零输出数组','Zero output',0,[]);
 if(kind==='adjoint'){const Rraw_b=R_b.slice();raw.constraints.forEach(i=>Rraw_b[i]=0);if(m.node)step('反传约束并屏蔽旧平衡行','Reverse constraints and mask old balance rows',0,raw.constraints.map(i=>['R_b['+i+']',R_b[i]]),L.constraintVjp(R_b).map((x,i)=>[i,x]));for(let f=m.m;f>=0;f--){const H_b=[0,1,2].map(k=>(f>0?Rraw_b[3*(f-1)+k]:0)-(f<m.m?Rraw_b[3*f+k]:0)),back=L.faceVjp(f,H_b);step('形成净面通量种子','Form net-face seeds',f,H_b.map((v,k)=>['H_b['+k+']',v]));if(f>0&&f<m.m){const z=raw.local[f],tau_b=-H_b[1]-z.u*H_b[2],q_b=H_b[2],u_b=-z.tau*H_b[2];step('反传黏性功与热通量','Reverse viscous work and heat flux',f,[['u_b (viscous)',u_b],['tau_b',tau_b],['q_b',q_b]]);step('反传应力、梯度与输运系数','Reverse stress, gradients and transport',f,[['mu_b',4/3*z.ux*tau_b],['ux_b',4/3*z.mu*tau_b],['kappa_b',-z.tx*q_b],['tx_b',-z.kappa*q_b]]);}step('反传通量、物性、梯度和重构，累加状态','Reverse flux, transport, gradients and reconstruction; accumulate state',f,back.map((v,i)=>['U_b['+i+'] +=',v]).filter(a=>Math.abs(a[1])>1e-15),back.map((v,i)=>[i,v]));}}
 else{for(let f=0;f<=m.m;f++){const z=kind==='tangent'?tan.local[f]:raw.local[f];step('收集并重构左右状态','Gather and reconstruct traces',f,[...z.ul.map((v,k)=>[nm('Uminus['+k+']'),v]),...z.ur.map((v,k)=>[nm('Uplus['+k+']'),v])]);step('计算梯度和输运系数','Evaluate gradients and transport',f,['ux','tx','mu','kappa'].map(k=>[nm(k),z[k]]));step('计算应力、热流及能量黏性通量','Evaluate stress, heat flux and viscous work',f,[[nm('tau'),z.tau],[nm('q'),z.heat],...z.fv.map((v,k)=>[nm('Fv['+k+']'),v])]);step('计算对流与净通量','Evaluate convective and net flux',f,[...z.fc.map((v,k)=>[nm('Fc['+k+']'),v]),...z.H.map((v,k)=>[nm('H['+k+']'),v])]);const updates=[];z.H.forEach((v,k)=>{if(f>0)updates.push([3*(f-1)+k,v]);if(f<m.m)updates.push([3*f+k,-v]);});step('向左右控制体分发','Scatter to neighbouring volumes',f,[],updates);}
 const src=kind==='tangent'?tan.sources:raw.sources;src.forEach((ss,i)=>step('减去体积源项','Subtract integrated source',i,ss.map((v,k)=>[nm('S['+k+']'),v]),ss.map((v,k)=>[3*i+k,-v])));if(m.node)step('替换边界约束行','Replace boundary rows',0,[],raw.constraints.map(i=>[i,kind==='tangent'?tan.R[i]:raw.R[i],true]));}
 const expected=kind==='adjoint'?L.vjp(R_b).U:kind==='tangent'?tan.R:raw.R;return {steps,expected,error:N.norm(out.map((v,i)=>v-expected[i]))};}
const api={reference,accuracy,trace};if(typeof module!=='undefined')module.exports=api;else root.NSStudy=api;
})(typeof window!=='undefined'?window:globalThis);