本页的方程是对流方程的格点页面上的方程加一项常系数扩散 −ν∇2u（ν=0.5），算例是同一张环形（O 型）混合网格。有了扩散，内圈物面才能施加 Dirichlet 壁面条件 u=uw=0；本页把它强施加：壁面节点的方程在边循环之后被整行换成 ui−uw。扩散通量要用节点梯度，这里用最小二乘来算，所以残差成了两遍循环：先算梯度，再算通量。目标函数是流进物面的总通量。对照阅读：同一方程的格心页面。
The equation is that of the node-centred advection page with a constant-coefficient diffusion term −ν∇2u (ν=0.5) added, on the same annular (O-type) hybrid mesh. With diffusion the inner body surface can carry a Dirichlet wall condition, u=uw=0, and this page imposes it strongly: after the edge loop, each wall node's equation is replaced as a whole by ui−uw. The diffusive flux needs nodal gradients, taken here by least squares, so the residual becomes two loops: gradients first, then fluxes. The objective is the total flux into the body. Read alongside: the cell-centred page for the same equation.
把 u(x,y,t) 看作一个被输运的标量。任取一块区域，其中 u 的总量只因穿过边界的通量和区域内的源项 s 而改变。通量有两部分：对流通量 12u2β，方向由常向量 β 给出，与对流页面第 1 节的方程相同；扩散通量 −ν∇u，从 u 高处指向低处，大小与梯度成正比，ν>0 是扩散系数。写成微分形式，就是
Think of u(x,y,t) as a transported scalar. The amount of u in any region changes only through the flux across its boundary and the source s inside it. The flux has two parts: the convective flux 12u2β, its direction set by the constant vector β, as in the equation of Section 1 of the advection page; and the diffusive flux −ν∇u, which points from high to low u with a size proportional to the gradient, ν>0 being the diffusion coefficient. In differential form,
按乘积法则（β 是常向量），∇·(12u2β)=uβ·∇u：u 沿 β 以正比于自身的速度被输运，与一维 Burgers 方程的 uux 同类。扩散项是 −∇·(ν∇u)=−ν∇2u。本页只研究定常解，基准问题也没有源项。令 ∂u/∂t=0、s=0，得到要离散的方程：
By the product rule (β is constant), ∇·(12u2β)=uβ·∇u: u is carried along β at a speed proportional to its own value, of the same kind as the uux of the one-dimensional Burgers equation. The diffusive term is −∇·(ν∇u)=−ν∇2u. This page studies steady solutions only, and the baseline problem has no source. Setting ∂u/∂t=0 and s=0 gives the equation to discretise:
第一项就是对流页面的非线性守恒律，对流通量为 12u2β；第二项是常系数扩散，扩散通量为 −ν∇u。求解域 Ω 是物体与远场之间的环形区域。源项并没有消失：每个控制体一个源项强度，作为设计变量留在离散方程里，第 11 节对它求导，第 15 节用它构造精确解。
The first term is the nonlinear conservation law of the advection page, with convective flux 12u2β; the second is constant-coefficient diffusion, with diffusive flux −ν∇u. The domain Ω is the annulus between the body and the far field. The source has not gone away: one source strength per control volume stays in the discrete equations as a design variable, Section 11 differentiates with respect to it, and Section 15 uses it to manufacture an exact solution.
图 1.1 把本页的问题画在一起：物体（琥珀色）与远场之间的求解域、两类边界，以及本页程序在这张网格上算出的解 u*（颜色越浓，u 越大）。u* 在物面上是 0，在迎着来流的远场附近最大；下面几小节逐一说明这张图。
Figure 1.1 puts the problem of this page in one picture: the domain between the body (amber) and the far field, its two kinds of boundary, and the solution u* that this page's code computes on this mesh (the stronger the colour, the larger u). u* is 0 on the body and largest near the part of the far field that faces the oncoming flow; the subsections below take the picture apart.
没有扩散时（对流页面），方程只含一阶导数，是双曲型的：通量对 u 的导数 uβ 是特征速度，在 u>0 的地方它指向 β 的方向，信息只沿这个方向传播，边界上只能在特征进入区域的地方给数据。沿着一个封闭的物体，β·n 必然有正有负：本网格物体的 8 条边中，5 条迎着来流（β 从流体穿进物体），3 条背着来流（图 1.1）。所以纯对流时不能在整条物面上给定 u：在迎着来流的那 5 条边上，u 已经由上游的数据决定了。
Without diffusion (the advection page) the equation has first derivatives only and is hyperbolic: the derivative of the flux with respect to u, uβ, is the characteristic speed, which points along β wherever u>0; information travels only in that direction, and boundary data can be given only where the characteristics enter the domain. Along a closed body β·n must take both signs: of the 8 sides of this mesh's body, 5 face the flow (β passes from the fluid into the body) and 3 face away from it (Figure 1.1). So pure advection cannot take a prescribed u along the whole body: on the 5 sides that face the flow, u is already set by the data upstream.
加上扩散项，方程成为二阶的。定常方程的最高阶项是 −ν∇2u，它是椭圆型的：一点的值受四周各个方向的影响，上游和下游都有，边界上每一处都要、也只要一个条件。本页的两类边界因此是：
With the diffusion term the equation becomes second order. The highest-order term of the steady equation is −ν∇2u, which makes it elliptic: the value at a point depends on its surroundings in every direction, upstream and downstream alike, and every part of the boundary takes exactly one condition. So the two kinds of boundary of this page are:
对流与扩散孰强孰弱，由 Péclet 数衡量，它是对流输运与扩散输运之比：Pe=|β|u∞L/ν。取物体的直径（两个壁面节点间的最大距离）L=1.232，|β|=1.25≈1.118，得 Pe≈2.75：对流强一些，但两者同一量级。在来流速度下，扩散与对流相平衡的长度是 ν/(|β|u∞)≈0.45，是物面到远场的距离（1.0 到 1.6）的 28% 到 45%。所以 u 从壁面的 0 升向来流值的那一层并不薄，它占了求解域的很大一部分（图 1.1）。在这张粗网格上，对流通量里的耗散也像一个与 ν 同量级的扩散系数（第 3 节末尾的注），这是 ν 取 0.5 而不更小的原因。
Which of the two is stronger is measured by the Péclet number, the ratio of convective to diffusive transport: Pe=|β|u∞L/ν. With the body's diameter (the largest distance between two wall nodes), L=1.232, and |β|=1.25≈1.118, Pe≈2.75: convection is the stronger, but the two are of the same order. At the free-stream speed, the length on which diffusion balances convection is ν/(|β|u∞)≈0.45, 28% to 45% of the distance from the body to the far field (1.0 to 1.6). So the layer in which u climbs from 0 at the wall towards the free-stream value is not thin: it fills a large part of the domain (Figure 1.1). On this coarse mesh the dissipation in the convective flux also acts like a diffusivity of the same size as ν (the note at the end of Section 3), which is why ν is 0.5 and not smaller.
图 1.2 把基准解上的两种通量逐点画出来。按每一圈的总和看，贴壁的一圈里只有扩散通量（壁面节点上 u=0），最外一圈里对流通量占上风；壁面节点上的扩散通量都指向物体。
Figure 1.2 draws the two fluxes of the baseline solution point by point. Summed ring by ring, the ring on the wall carries diffusion only (u=0 at the wall nodes) and convection wins in the outermost ring; at the wall nodes the diffusive flux points into the body.
物体是一个不规则的八边形，远场是半径 1.9 的圆的内接八边形，两者之间就是求解域，与对流页面相同。数据是 β=(1,0.5)、ν=0.5、来流值 u∞=1 与壁面值 uw=0。基准解 u* 在 8 个壁面节点上正好等于 uw=0（强施加，第 4 节），其余 16 个节点上介于 0.151 与 0.808 之间：最小的在背着来流一侧的中圈节点 n10，最大的在迎着来流的远场节点 n22（图 1.1；由第 6 节的 Jacobi 迭代算出）。u* 从物面向外增大，但处处都没有达到 u∞=1，因为远场条件不是 u=u∞（第 4 节）。
The body is an irregular octagon and the far field the octagon inscribed in the circle of radius 1.9; the domain lies between them, as on the advection page. The data are β=(1,0.5), ν=0.5, the free-stream value u∞=1 and the wall value uw=0. The baseline solution u* equals uw=0 exactly on the 8 wall nodes (imposed strongly, Section 4) and lies between 0.151 and 0.808 on the other 16: the smallest on the middle-ring node n10 on the side facing away from the flow, the largest on the far-field node n22 facing the oncoming flow (Figure 1.1; computed by the Jacobi iteration of Section 6). u* grows away from the body but nowhere reaches u∞=1, because the far-field condition is not u=u∞ (Section 4).
第 7 节的目标函数是流进物体的总通量 J，它从被替换的壁面行里找回来。在基准解上 J=2.748718：从远场的 8 个入流半面流入 2.907501，从 8 个出流半面流出 0.158784，两者之差全部流进物体。壁面上 u=0，对流通量 12u2β·n 为零，所以流进物体的全是扩散通量。
The objective of Section 7 is the total flux into the body, J, recovered from the replaced wall rows. At the baseline solution J=2.748718: 2.907501 flows in through the 8 inflow half-faces of the far field and 0.158784 out through the 8 outflow half-faces, and the difference all flows into the body. On the wall u=0, so the convective flux 12u2β·n vanishes there, and all the flux into the body is diffusive.
设计输入 D=(u∞,si,uw,i,X) 有四部分：来流值 u∞；每个控制体 Vi 里的源项强度 si（基准值为 0）；每个壁面节点上的壁面值 uw,i（基准值为 0）；以及 24 个网格节点的坐标 X（网格见第 2 节）。合起来是 1 + 24 + 8 + 48 = 81 个分量。按伴随优化文献的通常写法，设计变量记作 D。第 10、11 节求 J 对 u∞ 的导数，第 11 节求它对 si 与 uw,i 的导数，第 12 节求它对 X 的导数。
The design input D=(u∞,si,uw,i,X) has four parts: the free-stream value u∞; a source strength si in each control volume Vi (0 in the baseline); a wall value uw,i at each wall node (0 in the baseline); and the coordinates X of the 24 mesh nodes (the mesh is Section 2's). Together that is 1 + 24 + 8 + 48 = 81 components. As is usual in the adjoint-optimization literature, the design variables are written D. Sections 10 and 11 take the derivative of J with respect to u∞, Section 11 those with respect to si and uw,i, and Section 12 those with respect to X.
在流体力学里，对流扩散方程描述一个量同时被流动携带（对流）、又因分子运动而扩散的过程：温度（热传导）与浓度（物质扩散）都满足这类方程，动量方程里对流项与黏性项之间也是同样的平衡。本页的对流通量 12u2β 是非线性的，输运速度 uβ 由 u 自己决定，所以它更接近动量方程：一维时它就是带黏性的 Burgers 方程。物面上的 u=uw 相当于黏性流动的无滑移条件；流进物体的总通量则相当于壁面热流或壁面摩擦，都是要用壁面上的梯度才能算出的量。
In fluid dynamics the advection–diffusion equation describes a quantity that is carried by the flow (advection) and spreads by molecular diffusion at the same time: temperature (heat conduction) and concentration (mass diffusion) obey equations of this kind, and in the momentum equation the convection and viscous terms strike the same balance. The convective flux here, 12u2β, is nonlinear, its transport speed uβ set by u itself, which brings it closer to the momentum equation: in one dimension it is the viscous Burgers equation. The wall condition u=uw plays the part of the no-slip condition of a viscous flow, and the total flux into the body that of the wall heat flux or the skin friction, quantities that need the gradient at the wall.
在本课程里，本页是从纯对流走向 Navier–Stokes 方程黏性项的一步。对流页面的面通量只要两侧的值；本页的扩散通量要面上的梯度，于是先用一遍循环算出最小二乘梯度，再算通量（第 5 节的两遍循环），精确 Jacobian 伸到邻居的邻居（第 6 节），物面可以施加 Dirichlet 条件（第 4 节）。二维层流 N–S 与 RANS–SA 页面的黏性应力和热通量也要面上的梯度，RANS–SA 页面的面梯度正是按本页第 3 节的办法构造的。未知量仍是每点一个标量，这些新东西都能在 24 个未知量上一一看清。
In this course the page is the step from pure advection towards the viscous terms of the Navier–Stokes equations. The advection page's face flux needs only the values on either side; the diffusive flux here needs the gradient on the face, so a first loop computes least-squares gradients before a second computes the fluxes (the two loops of Section 5), the exact Jacobian reaches the neighbours of neighbours (Section 6), and the body can take a Dirichlet condition (Section 4). The viscous stresses and heat fluxes of the 2-D laminar N–S and RANS–SA pages need gradients on the faces too, and the RANS–SA page builds its face gradient just as Section 3 of this page does. The unknown is still one scalar per point, so all of this can be followed on 24 unknowns.
图 1.3 把本页放回网站的路线里：上排是标量页面，从一维 Burgers 到二维对流，再到本页；下排是方程组页面。点击任一方框，进入那一页的第 1 节。
Figure 1.3 puts this page on the route through the site: the top row holds the scalar pages, from 1-D Burgers to 2-D advection and on to this page; the bottom row the systems. Click any box to go to that page's Section 1.
本页要回答的问题是：流进物体的总通量 J 怎样随设计输入 D 的 81 个分量中的每一个变化，也就是梯度 dJ/dD。前向（tangent）方法每个分量解一次线性方程，共 81 次；伴随方法只解一次转置方程，再做一次乘积，就得到全部 81 个分量。页首的全页路线图把这两条路线并排画出，方框下的节号指向讲它们的各节。
The question this page answers: how does the total flux into the body, J, change with each of the 81 components of the design input D, that is, what is the gradient dJ/dD? The tangent (forward) method solves one linear system per component, 81 in all; the adjoint method solves one transposed system and gets all 81 components with one more product. The page roadmap at the top draws the two routes side by side, and the section numbers under its boxes lead to the sections that treat them.
要回答它，先得把方程变成有限个数。第 2 节选网格、定下未知量放在哪里，第 3 节在每个对偶面上写出对流通量与扩散通量。
To answer it, the equation must first become finitely many numbers. Section 2 chooses the mesh and where the unknowns sit, and Section 3 writes the convective and the diffusive flux on each dual face.
第 1 节的方程在求解域的每一点都成立，而计算机只能存有限个数。本节先把求解域分成网格：它与对流页面是同一张混合网格，未知量放在节点上，控制体是中位对偶体；扩散还要用到两样新的几何量和每个节点的梯度模板。
The equation of Section 1 holds at every point of the domain, but a computer stores finitely many numbers. This section divides the domain into a mesh: the same hybrid mesh as the advection page, with the unknowns on the nodes and median duals as control volumes; diffusion also needs two new geometric quantities and the gradient stencil of each node.
网格与对流页面第 2 节完全相同：3 圈节点、每圈 8 个，共 24 个节点；贴壁一层是 8 个四边形，外层是 16 个三角形。未知量放在节点上，控制体是中位对偶体。它的构造方法（在每个单元内把边中点与单元中心相连）在对流页面里有逐步演示，这里不再重复。网格有 48 条边，对应 48 个内部对偶面；另有 32 个边界半面（16 个在壁面上，16 个在远场上），所以边循环一共要走 80 个面。
The mesh is exactly the one of Section 2 of the advection page: 3 rings of 8 nodes, 24 nodes in all, with a layer of 8 quadrilaterals against the wall and 16 triangles outside it. The unknowns sit on the nodes and the control volumes are the median duals. How those are built (inside every cell, join the edge midpoints to the cell centre) is shown step by step on the advection page and is not repeated here. The 48 edges give 48 interior dual faces, and there are 32 boundary half-faces besides (16 on the wall, 16 on the far field), so the edge loop visits 80 faces.
变的是节点的角色（图 2.1）。内圈 8 个节点落在物面上，而物面现在是 Dirichlet 壁面，这些节点的值是给定的 uw。它们仍然算未知量，在 24 个里占 8 个位置，只是它们的方程要被换掉（第 4 节）。外圈 8 个节点在远场上，其余 8 个是内部节点。
What changes is the role of the nodes (Figure 2.1). The 8 inner-ring nodes lie on the body, which is now a Dirichlet wall, so their values are prescribed: uw. They still count as unknowns and take 8 of the 24 slots; only their equations will be swapped (Section 4). The 8 outer-ring nodes lie on the far field, and the remaining 8 are interior nodes.
对流页面的边循环在每个面上只需要两样几何量：法向 n 与它的长度 |n|。扩散项还要两样。第一样是边向量 ξ=xj−xi，它连接对偶面两侧的两个未知量，两点差 uj−ui 就是沿它量出来的。第二样是每个节点的梯度模板，即最小二乘梯度要用哪些邻居的值。本页取最简单的选择：边邻居，也就是和该节点共用一条边的节点。壁面节点有 3 个边邻居，中圈节点有 5 个，外圈节点有 4 个。
The edge loop of the advection page needed two geometric quantities per face: the normal n and its length |n|. Diffusion needs two more. The first is the edge vector ξ=xj−xi, which joins the two unknowns on either side of the dual face; the two-point difference uj−ui is measured along it. The second is each node's gradient stencil, the neighbours whose values the least-squares gradient uses. This page takes the simplest choice, the edge neighbours: the nodes that share an edge with the node. A wall node has 3 of them, a middle-ring node 5, an outer-ring node 4.
在 图 2.2 里点选任意节点，可以看到它的控制体、对偶面和梯度模板（虚线箭头就是各个 ξs）；再点选一个对偶面，可以看到它的 n、ξ，以及下一节要用的分解 n=pξ+k。
In Figure 2.2, pick any node to see its control volume, its dual faces and its gradient stencil (the dashed arrows are the ξs); then pick a dual face to see its n and ξ, and the split n=pξ+k used in the next section.
网格、未知量的位置和梯度模板都定下来了。第 3 节把第 1 节的方程在每个控制体上积分，写出每个对偶面上的对流通量与扩散通量，以及扩散通量要用的最小二乘梯度。
The mesh, the place of the unknowns and the gradient stencils are now fixed. Section 3 integrates the equation of Section 1 over each control volume and writes the convective and the diffusive flux on every dual face, with the least-squares gradients the diffusive flux needs.
本节写出每个对偶面上的通量：先把方程化成各个对偶面上通量之和（有限体积离散），再写对流通量（与对流页面相同）和扩散通量；扩散通量要面上的法向导数，因此还要最小二乘梯度。
This section writes the flux on every dual face: first the equation becomes a sum of dual-face fluxes (the finite-volume discretisation), then come the convective flux (as on the advection page) and the diffusive flux, whose normal derivative on the face needs least-squares gradients.
把第 1 节的方程在节点 i 的控制体 Vi 上积分，再用散度定理把体积分化成边界积分：
Integrate the equation of Section 1 over the control volume Vi of node i and turn the volume integral into a boundary integral with the divergence theorem:
求和取遍节点 i 的所有对偶面 e。每个面上的通量只算一次，然后以相反符号分给两侧：i 一侧取 +，j 一侧取 −，所以离散格式是守恒的。ha 是对流通量（下标 a 取自 advection），hv 是扩散通量。这里的 n 是单位法向在整个对偶面上的积分（由 i 指向 j，边界半面上指向域外），所以 vn、ha、hv 都已含面的大小；上面积分里的 n 则是单位法向。有源项时，Ri 再减去控制体 Vi 里的源项强度：Ri←Ri−si；基准问题里 si=0，它作为设计变量留在方程里（第 11 节）。
The sum runs over the dual faces e of node i. Each face's flux is computed once and handed to its two sides with opposite signs, + for side i and − for side j, so the scheme is conservative. ha is the convective flux (the subscript is from advection) and hv the diffusive one. Here n is the unit normal integrated over the whole dual face (from i to j, out of the domain on a boundary half-face), so vn, ha and hv already carry the face's size; in the integral above n is the unit normal. With a source, Ri also subtracts the source strength si of Vi, Ri←Ri−si; si=0 in the baseline, and it stays in the equations as a design variable (Section 11).
中心平均加常数耗散，外加远场上由 vn 的符号冻结的入流、出流开关。它的推导在对流页面第 3、4 节，这里原样沿用。
A central average plus constant dissipation, with the inflow/outflow switch on the far field frozen by the sign of vn. It is derived in Sections 3 and 4 of the advection page and is used here unchanged.
扩散通量 hv=−ν∇u·n 要的是面上的法向导数。面两侧只有两个节点值 ui、uj，而二维梯度有两个分量，只靠这两个值不够；除非对偶面恰好垂直于边（n 与 ξ 平行），那时两点差就是法向导数。一般网格上 n 与 ξ 并不平行（图 3.1），做法是把 n 分成沿 ξ 的部分和垂直于 ξ 的部分：
The diffusive flux hv=−ν∇u·n needs the normal derivative on the face. Only two nodal values sit either side of it, ui and uj, while a two-dimensional gradient has two components, so those two values are not enough, unless the dual face happens to be perpendicular to the edge (n parallel to ξ), in which case the two-point difference is the normal derivative. On a general mesh n and ξ are not parallel (Figure 3.1), and the way out is to split n into a part along ξ and a part perpendicular to it:
沿 ξ 的部分用两点差 ∇u·ξ≈uj−ui；垂直的部分 k 用两端节点梯度的平均 Ge=12(Gi+Gj)：
The part along ξ uses the two-point difference ∇u·ξ≈uj−ui; the perpendicular part k uses the average of the two nodal gradients, Ge=12(Gi+Gj):
它也常写成先修正平均梯度、再投影到 n 上：Ge·n−(Ge·ξ−Δu)p，展开就是上式。两点部分让 i 与 j 直接耦合，是格式稳定的来源；梯度部分补上 n 偏离 ξ 的那一块。两部分对线性场都精确，所以整个 hv 对线性场精确。
It is often written as a correction of the averaged gradient before projecting it on n, Ge·n−(Ge·ξ−Δu)p; expanded, that is the line above. The two-point part couples i and j directly and is where the scheme's stability comes from; the gradient part makes up for n leaning away from ξ. Both parts are exact on linear fields, so the whole hv is.
节点 i 的梯度 Gi 取为让线性模型 ui+G·ξs 最贴近模板中各邻居值 us 的那个 G（ξs=xs−xi）：
The gradient Gi of node i is the G for which the linear model ui+G·ξs comes closest to the neighbour values us of the stencil (ξs=xs−xi):
权重 wis 只依赖几何，每张网格算一次。因为 ΣswisξsT=Mi−1Mi=I，线性场的梯度被精确复现。本页不加距离权（非加权最小二乘）。梯度对 u 是线性的，这一点在第 8 节会用到。在 图 3.3 里可以选一个节点，看它的模板、矩阵 Mi、各个权重，以及对三种场取出的梯度。
The weights wis depend on the geometry only and are computed once per mesh. Because ΣswisξsT=Mi−1Mi=I, the gradient of a linear field comes out exactly. This page uses no distance weighting (unweighted least squares). The gradient is linear in u, which Section 8 will use. In Figure 3.3, pick a node to see its stencil, its matrix Mi, each weight, and the gradient taken from three fields.
上面的公式是一个最小二乘问题的解。下面从问题出发，一步一步推到公式，再推到第 6 节的导数和第 9 节的转置。
The formula above is the solution of a least-squares problem. The steps below start from the problem, arrive at the formula, and go on to the derivative that Section 6 needs and the transpose that Section 9 needs.
第一步，写出要最小化的量。在节点 i 的每个边邻居 s 上，线性模型 ui+G·ξs 与那里的值 us 差 us−ui−G·ξs；最小二乘取使这些差的平方和最小的 G：
Step 1, the quantity to minimise. At each edge neighbour s of node i the linear model ui+G·ξs misses the value us there by us−ui−G·ξs; least squares takes the G that makes the sum of their squares smallest:
第二步，令它对 G 的导数为零。平方和对 G 求导得 −2Σsξs(us−ui−G·ξs)；再用 ξs(G·ξs)=(ξsξsT)G，就得到 2×2 的正规方程：
Step 2, set its derivative with respect to G to zero. The derivative of the sum of squares is −2Σsξs(us−ui−G·ξs); with ξs(G·ξs)=(ξsξsT)G this gives the 2×2 normal equations:
第三步，解出 Gi。左边的 2×2 矩阵 Mi（正文的记号）只依赖几何。由 Cauchy–Schwarz 不等式，它的行列式不小于零，只有当所有 ξs 都平行时才等于零。两边乘它的逆，就是正文的公式：
Step 3, solve for Gi. The 2×2 matrix Mi on the left (the text's symbol) depends on the geometry only. By the Cauchy–Schwarz inequality its determinant is never negative, and it is zero only when all the ξs are parallel. Multiplying by its inverse gives the formula of the text:
2×2 矩阵的逆可以直接写出来，代码 lsq_weights 正是这样算的，每张网格只算一次：
The inverse of a 2×2 matrix can be written out, and that is exactly how the code, lsq_weights, computes it, once per mesh:
第四步，线性场被精确复现。设 u 是梯度为常向量 ∇u 的线性场，则 us−ui=ξs·∇u=ξsT∇u，代入第三步：
Step 4, linear fields come out exactly. Let u be a linear field with a constant gradient ∇u; then us−ui=ξs·∇u=ξsT∇u, and step 3 gives
这就是 ΣswisξsT=Mi−1Mi=I 的意思。
That is what ΣswisξsT=Mi−1Mi=I says.
第五步，对状态线性化。Gi 对 u 是线性的，偏导数逐项读出（m 不在模板里时 wim=0）：
Step 5, the linearisation in the state. Gi is linear in u, and its partial derivatives can be read off term by term (wim=0 when m is not in the stencil):
这就是第 6 节 ∂Gi/∂us=wis、∂Gi/∂ui=−Σswis 的来历，也是代码 gradRows() 存下的量。
This is where the ∂Gi/∂us=wis and ∂Gi/∂ui=−Σswis of Section 6 come from, and what the code stores in gradRows().
第六步，转置：伴随怎样把第一遍倒过来。（记号：加点的 ẋ 是前向（tangent）模式的导数，加横线的 x̄ 是反向（adjoint）模式累加的伴随，第 8 节的“记号”一段有完整说明。）梯度的伴随 Ḡi 与梯度的切线配成下面的对偶恒等式；按状态的切线 u̇ 读出系数，就是状态的伴随 ū 分到的份额：每个模板点得到一份、节点自己减去一份：
Step 6, the transpose: how the adjoint runs pass 1 backwards. (Notation: a dot, ẋ, is a forward (tangent) derivative and a bar, x̄, the adjoint accumulated in reverse mode; the Notation paragraph of Section 8 says more.) Pair the gradient's adjoint Ḡi with the gradient's tangent in the duality identity below; read off by the state tangents u̇, it gives the state adjoint ū its shares: every stencil point receives one and the node itself gives one up:
每条边 (i,j) 在两个模板里各出现一次（j 在 i 的模板里，i 在 j 的模板里）。按边组织这个循环，就是 residual_b 第 12 至 15 行的 γi=wij·Ḡi、γj=wji·Ḡj（代码把 Ḡ 存在数组 G_b 里，ū 存在 y 里）。
Every edge (i,j) appears in two stencils (j in that of i, i in that of j). Organised by edges, this loop is lines 12 to 15 of residual_b, γi=wij·Ḡi and γj=wji·Ḡj (the code keeps Ḡ in the array G_b and ū in y).
核对：本页 24 个节点上，按第三步的闭式算出的权重与代码的 wis 最多差 4.4×10−16；ΣswisξsT−I 的最大分量是 4.4×10−16；线性场 u=0.3+0.7x−0.4y 的最小二乘梯度与 (0.7,−0.4) 最多差 3.3×10−16；gradRows() 给出的 ∂G/∂u 与对 grads() 的中心差分最多差 5.7×10−14（本页代码在构建时计算）。
Check: over the 24 nodes of this page, the weights from the closed form of step 3 differ from the code's wis by at most 4.4×10−16; the largest entry of ΣswisξsT−I is 4.4×10−16; for the linear field u=0.3+0.7x−0.4y the least-squares gradient misses (0.7,−0.4) by at most 3.3×10−16; and the ∂G/∂u of gradRows() differs from central differences of grads() by at most 5.7×10−14 (computed by this page's code at build time).
图 3.2 画出三种节点的模板向量 ξs 与权重 wis：权重一般不与 ξs 同向，正是 2×2 矩阵的逆把它们调成对线性场精确的组合。
Figure 3.2 draws the stencil vectors ξs and the weights wis of three kinds of node: a weight does not in general point along its ξs, and it is the inverse of the 2×2 matrix that turns them into the combination that is exact on linear fields.
选线性场时，任何节点的误差都在舍入量级；选二次场时，误差随节点间距变化，这是一阶近似的正常表现。选收敛解 u* 时看到的，就是第二遍（通量循环）真正使用的梯度。
With the linear field the error is at round-off level on every node; with the quadratic one it varies with the node spacing, as a first-order approximation should. With the converged solution u* you see exactly the gradients that pass 2, the flux loop, will use.
图 3.4 把一个内部对偶面上的量全部列出来，状态取收敛解 u*。表里同时给出 ∇u·n 的三种估计：只用平均梯度 Ge·n、只用两点差 pΔu，以及修正后的 pΔu+Ge·k。面越歪（|k|/|n| 越大），只用两点差的估计离修正值越远，因为两者之差正是 Ge·k；只用平均梯度的估计差 p(Ge·ξ−Δu)，它与面歪不歪无关，在几乎正交的面上也可能很大。
Figure 3.4 lists everything on one interior dual face, at the converged solution u*. The table also gives three estimates of ∇u·n: the averaged gradient alone, Ge·n; the two-point difference alone, pΔu; and the corrected pΔu+Ge·k. The more skewed the face (the larger |k|/|n|), the further the two-point difference alone lies from the corrected value, since the two differ by exactly Ge·k; the averaged gradient alone misses by p(Ge·ξ−Δu), which does not depend on the skewness and can be large even on a nearly orthogonal face.
第一步，p 由“k 与 ξ 垂直”唯一确定：
Step 1, p is fixed by requiring k to be perpendicular to ξ:
第二步，对线性场精确。设 u 是梯度为常向量 ∇u 的线性场。ξ 连接 uL、uR 所在的两点，所以 uR−uL=ξ·∇u；由上一段第四步，两端的最小二乘梯度都等于 ∇u，所以平均梯度 Ge=∇u。于是
Step 2, exactness on linear fields. Let u be a linear field with a constant gradient ∇u. ξ joins the two points where uL and uR live, so uR−uL=ξ·∇u; by step 4 above, the least-squares gradients on both sides equal ∇u, so the averaged gradient Ge=∇u. Hence
两点部分和梯度部分各管 n 的一段，加起来正好是 n。
The two-point part and the gradient part each take one piece of n, and together they take all of it.
第三步，对四个输入线性化：两侧的状态 uL、uR 与两侧的梯度 GL、GR。对流部分 ha=14vn(uL2+uR2)−ε|n|(uR−uL) 只含两个状态；扩散部分对四个输入都是线性的：
Step 3, linearise with respect to the four inputs: the two states uL, uR and the two gradients GL, GR. The convective part ha=14vn(uL2+uR2)−ε|n|(uR−uL) holds the two states only; the diffusive part is linear in all four:
第四步，这些偏导数就是 residual_d 第 12 行 ḣ=cLu̇L+cRu̇R−ν[p(u̇R−u̇L)+Ġe·k] 里的系数，也是 residual_b 反过来乘的系数。再经第一遍链式求导（上一段第五步），就得到第 6 节的公式（[m=L] 在 m=L 时为 1，否则为 0）：
Step 4, these partial derivatives are the coefficients on line 12 of residual_d, ḣ=cLu̇L+cRu̇R−ν[p(u̇R−u̇L)+Ġe·k], and the ones residual_b multiplies by on the way back. Chained through pass 1 (step 5 above) they give the formula of Section 6 ([m=L] is 1 if m=L and 0 otherwise):
核对：在本页 48 个内部对偶面上，对 faceFlux() 取中心差分（步长 10−6），与上面的偏导数最多差 8.3×10−11；取线性场 u=0.3+0.7x−0.4y，修正后的扩散通量与 −νn·∇u 最多差 1.1×10−16（本页代码在构建时计算）。
Check: on the 48 interior dual faces of this page, central differences of faceFlux() (step 10−6) differ from the partial derivatives above by at most 8.3×10−11; for the linear field u=0.3+0.7x−0.4y, the corrected diffusive flux differs from −νn·∇u by at most 1.1×10−16 (computed by this page's code at build time).
图 3.5 把一个面的两遍循环画成一张计算图，边上正是这些偏导数。前向沿箭头走；伴随把同一张图倒过来走，每条边乘同一个数：这就是第 9 节 residual_b 的全部内容。
Figure 3.5 draws the two passes of one face as a computational graph whose edges carry exactly these partial derivatives. The forward mode follows the arrows; the adjoint walks the same graph backwards, multiplying by the same number on every edge: that is all there is to residual_b in Section 9.
v、G_d、h_d、r_d 与 y、G_b、h_b、w 里。数值取自本页的 u*。Figure 3.5 · The computational graph of the flux h on the interior face e20 (L=n9, R=n18). Boxes are stored quantities, and each arrow carries a local partial derivative. Top: forward (residual_d): the states' tangent u̇ enters the gradients' tangents Ġ through the least-squares weights in pass 1, enters ḣ in pass 2, and is handed to the two rows ṙ with +1 and −1. Bottom: adjoint (residual_b): the same graph with the same weights and every arrow reversed: the two rows' adjoints first combine into h̄=r̄L−r̄R, which then travels back along each edge, multiplied by the same partial derivative; a quantity with several outgoing arrows collects several terms with += on the way back. In the code these live in the arrays v, G_d, h_d, r_d and y, G_b, h_b, w. Numbers at this page's u*.ν 为什么取 0.5。耗散项 −ε|n|(uR−uL) 与扩散的两点部分 −νp(uR−uL) 形式相同，它就像一个数值扩散系数 ε|n|/p。在这张只有 24 个节点的粗网格上，这个系数在各内部面上介于 0.24 与 1.58 之间，每加密一次减半。在收敛解上，ν=0.5 时各内部面耗散通量的绝对值之和是扩散通量的 1.16 倍，两者相当；若取 ν=0.1，这个比值是 5.6，物理扩散只是数值耗散上的一点修正。第 15 节的精度阶研究会再遇到它。
Why ν=0.5. The dissipation −ε|n|(uR−uL) has the same form as the two-point part of diffusion, −νp(uR−uL), so it acts like a numerical diffusivity ε|n|/p. On this 24-node coarse mesh that diffusivity lies between 0.24 and 1.58 over the interior faces, and it halves with every refinement. At the converged solution with ν=0.5, the absolute dissipative fluxes over the interior faces sum to 1.16 times the diffusive ones: the two are comparable. With ν=0.1 the ratio would be 5.6, and the physical diffusion only a small correction to the numerical one. The order study of Section 15 meets it again.
到这里，内部对偶面上的通量都有了。边界还差两处：壁面节点的值给定为 uw，它们的方程要在边循环之后被替换；远场半面外侧没有节点，外侧状态要由来流值或外推补上。第 4 节处理这两处。
By now every interior dual face has its flux. Two things at the boundary are still missing: the wall nodes have the given value uw, and their equations are replaced after the edge loop; a far-field half-face has no node beyond it, and its outer state comes from the free stream or from extrapolation. Section 4 handles both.
物面上 u=uw=0。先看壁面节点 n1 的控制体（图 4.1）：它有 3 个内部对偶面，分别通往边邻居 n2、n8、n9，另外还有两个壁面半面。如果照常写守恒方程，就是“所有这些面的通量之和为零”。壁面半面上的通量是物面吸收（或放出）的量，事先并不知道：它就是让 u1 保持为 uw 所需要的那个量。由于 uw=0，对流通量 12u2β·n 在物面上为零，穿过物面的只有扩散通量。
On the body u=uw=0. Look first at the control volume of wall node n1 (Figure 4.1): it has 3 interior dual faces, towards its edge neighbours n2, n8 and n9, and two wall half-faces. Written as usual, its conservation equation says that the fluxes through all these faces sum to zero. The flux through the wall half-faces is what the body absorbs (or gives off), and it is not known in advance: it is whatever it takes to hold u1 at uw. Because uw=0, the convective flux 12u2β·n vanishes on the body, and only diffusion crosses it.
Dirichlet 条件有两种施加方式。弱施加：把 uw 写进壁面半面上的通量公式，方程照常组装；格心页面只能这样做，因为那里物面上没有未知量。强施加：不写壁面节点的守恒方程，而是把这一行整个换成
A Dirichlet condition can be imposed in two ways. Weakly: write uw into the flux formula of the wall half-faces and assemble the equations as usual; the cell-centred page can only do this, since it has no unknown on the body. Strongly: do not write the wall nodes' conservation equations at all, and replace each of those rows as a whole by
本页用强施加。它的好处是 ui=uw 在收敛解上严格成立，而不只是在网格加密的极限下成立，而且那个未知的壁面通量根本用不到。代价有两个：壁面通量要事后找回来（第 7 节）；前向和伴随都要为这几行做专门的处理（第 8、9 节）。
This page imposes the wall strongly. The gain is that ui=uw holds exactly at the converged solution, not merely in the limit of refinement, and the unknown wall flux is never needed. There are two prices: the wall flux has to be recovered afterwards (Section 7), and the forward and the adjoint code both need special treatment of these rows (Sections 8 and 9).
实现很简单。边循环照常把所有内部对偶面与远场半面的通量累加进 r，壁面半面直接跳过：它们的通量只会进入壁面节点的那一行，而那一行马上就要被丢掉。循环结束后，对每个壁面节点令 ri←ui−uw。在 Jacobian 里，这一行变成单位行 êiT：对角元为 1，其余为 0。替换后的残差与 Jacobian 记为 Rbc 和 Abc；替换之前的记为 R 和 A，第 7 节还要用到它们。
The implementation is short. The edge loop accumulates the fluxes of all interior dual faces and far-field half-faces into r as usual and simply skips the wall half-faces: their fluxes would only go into the wall nodes' rows, which are about to be thrown away. After the loop, set ri←ui−uw for every wall node. In the Jacobian this row becomes the unit row êiT: 1 on the diagonal, 0 elsewhere. The replaced residual and Jacobian are written Rbc and Abc; the ones before the replacement are R and A, and Section 7 still needs them.
远场半面上只算对流通量，扩散通量取零。对流部分与对流页面相同：远场半面按 vn=β·n 的符号分成入流和出流（在初始几何上冻结），入流取右状态 uR=u∞=1，出流取 uR=uL。16 个远场半面中，8 个入流，8 个出流。图 4.2 把每个远场半面上的通量拆开：中心对流部分 14vn(uL2+uR2)，耗散部分 −ε|n|(uR−uL)，扩散部分 0。
The far-field half-faces carry the convective flux only; their diffusive flux is set to zero. The convective part is as on the advection page: the far-field half-faces are split into inflow and outflow by the sign of vn=β·n (frozen at the initial geometry), the right state being uR=u∞=1 on inflow and uR=uL on outflow. Of the 16 far-field half-faces, 8 are inflow and 8 outflow. Figure 4.2 takes the flux on every far-field half-face apart: the central convective part 14vn(uL2+uR2), the dissipative part −ε|n|(uR−uL), and the diffusive part, 0.
“扩散通量取零”并不等于连续问题里的 ∂nu=0。要看它对应哪个连续边界条件，得看远场半面上的总通量。把通量除以面长 |n|，记单位外法向 n̂=n/|n|，b=β·n̂；连续问题在边界上的总通量是 12bu2−ν∂nu，离散闭合规定的正是它：
"Zero diffusive flux" is not the continuous condition ∂nu=0. Which continuous boundary condition it stands for is decided by the total flux on the far-field half-faces. Divide the flux by the face length |n| and write n̂=n/|n| for the unit outward normal and b=β·n̂; the continuous problem's total flux through the boundary is 12bu2−ν∂nu, and that is exactly what the discrete closure prescribes:
出流一侧（b≥0）两边的对流项相消，剩下 ∂nu=0，是齐次 Neumann 条件。入流一侧（b<0）剩下的是一个非线性 Robin 条件：
On the outflow side (b≥0) the convective terms cancel and ∂nu=0 is left: a homogeneous Neumann condition. On the inflow side (b<0) what is left is a nonlinear Robin condition:
它既不是 u=u∞，也不是 ∂nu=0。在连续读法里，扩散通量的位置由两样东西填上：中心平均与边界状态自己的对流通量 12bu2 之差，加上耗散项。耗散系数在这里不随网格变小：它乘的是面长 |n|，除以面长之后留下 ε；随网格变小的只是内部面上的数值扩散（第 15 节）。
It is neither u=u∞ nor ∂nu=0. In the continuous reading the place of the diffusive flux is taken by two things: the difference between the central average and the boundary state's own convective flux 12bu2, and the dissipation. Nor does the dissipation coefficient shrink with the mesh here: it multiplies the face length |n|, and dividing by that length leaves ε; what shrinks with the mesh is only the numerical diffusion on the interior faces (Section 15).
本页的解正是这样。在本页网格上，两个远场半面都是入流的 3 个节点上，u 介于 0.723 与 0.808 之间，而不是 1。一致加密（第 15 节的网格族）到第 5 层，这样的节点上的 u 仍介于 0.68 与 0.97 之间；由最小二乘梯度量出的外法向导数 ∂nu 介于 0.07 与 0.54 之间，并不趋于零，而且与上式由 u 算出的值越来越近，两者之差在各个入流节点上的均方根逐层是 0.292、0.111、0.045、0.020、0.009、0.004。出流一侧的 ∂nu 则趋于零：它的最大绝对值从第 1 层起逐层是 0.084、0.074、0.054、0.034、0.020（图 4.3）。附录 A.1 第 3 小节正是从这个总通量条件出发，推导远场的伴随边界条件。
The page's solution does exactly this. On the page's mesh, at the 3 nodes whose two far-field half-faces are both inflow, u lies between 0.723 and 0.808, not at 1. Refined uniformly (the family of Section 15) to level 5, u at such nodes still lies between 0.68 and 0.97; the outward normal derivative ∂nu measured by the least-squares gradient lies between 0.07 and 0.54, does not go to zero, and approaches the value the condition above gives from u, the root mean square of the difference over the inflow nodes being 0.292, 0.111, 0.045, 0.020, 0.009 and 0.004 level by level. On the outflow side ∂nu does go to zero: its largest magnitude is 0.084, 0.074, 0.054, 0.034 and 0.020 from level 1 on (Figure 4.3). Appendix A.1, part 3, starts from exactly this total-flux condition to derive the adjoint's far-field boundary condition.
为什么对流页面没有这种壁面。纯对流方程是一阶的，边界上只能在特征进入区域的地方给数据；物面是闭合曲线，特征速度沿物面必然有正有负，在特征离开区域的那部分物面上再规定 u，就多给了一个条件（见对流页面第 4 节）。扩散项是二阶的，它让整条物面上都可以、也必须给一个条件，Dirichlet 壁面由此成为适定的。
Why the advection page has no such wall. A purely convective equation is first order: boundary data can be given only where characteristics enter the domain. The body is a closed curve, the characteristic speed takes both signs along it, and prescribing u where characteristics leave the domain is one condition too many (see Section 4 of the advection page). The diffusion term is second order: it lets, and requires, the whole body take one condition, and that is what makes the Dirichlet wall well posed.
本页每个节点存储 1 个未知量。以下用 表示把所有位置按代码顺序排成的列向量（标量问题每个位置只有一个分量）。它是一个有限维向量，不是连续函数。物理方程中的空间导数只用于说明模型；求离散伴随时，求导对象是本页实际装配的有限维函数。
Each node stores 1 unknown. Below, stacks all locations in code order (a scalar problem has one component per location). It is a finite-dimensional vector, not a continuous function. Spatial derivatives in the physical equations explain the model; discrete differentiation acts on the actual finite-dimensional functions assembled on this page.
是本页全部设计变量排成的向量：来流值 、各控制体的源项强度 、各壁面节点的壁面值 与全部节点坐标 。第 10 节以来流值 为例输入（末尾一块也逐类讲了其余的设计分量），第 11 节讲源项 与壁面值 ，第 12 节把 写成 ，例如 ；第 7 节的 就是这里的状态梯度。
stacks all design variables of this page: the free-stream value , the source strength of each control volume, the wall value of each wall node and all node coordinates . Section 10 takes the free-stream value as its example input (its last block goes through the other kinds of design component), Section 11 the sources and the wall values , and Section 12 writes as , e.g. ; the of Section 7 is the state gradient here.
这里 、 都是列向量，它们的转置才是偏导数行。求导采用自动微分（AD）的记号：程序变量 上加点的 是前向（tangent）模式沿给定种子方向的导数，加横线的 是反向（adjoint）模式累加到 上的权重，例如 与 。两者都是对本页有限维代码求导；变分符号 只用于连续伴随附录里连续函数的变分。
Both and are columns; their transposes are the derivative rows. Derivatives use the notation of automatic differentiation (AD): a dot on a program variable is its forward (tangent) derivative along the chosen seed direction, and a bar is the weight accumulated into in reverse (adjoint) mode, for example and . Both differentiate the finite-dimensional code of this page; the variation symbol is kept for variations of continuous functions in the continuous-adjoint appendix.
对流页面的残差是一个边循环：每个面收集两侧状态、算一维通量、分发到两侧。扩散通量要用面两端的节点梯度，而节点梯度又要用邻居的值，所以不能在同一遍循环里边算边用：必须先把所有节点的梯度算完，才能开始算通量。残差因此成了两遍循环，再加最后的替换：
The residual of the advection page is one edge loop: each face gathers the two states, computes the one-dimensional flux and scatters it to both sides. The diffusive flux needs the nodal gradients at both ends of the face, and a nodal gradient needs the neighbours' values, so the gradients cannot be computed and used in the same sweep: the gradients of all nodes must be finished before the first flux is computed. The residual therefore becomes two loops, plus the replacement at the end:
图 5.1 把这三步画成数据流，并且和第 8、9 节的前向、伴随放在一起：前向走同样的顺序，伴随把顺序整个倒过来。
Figure 5.1 draws the three steps as a data flow, next to the forward and adjoint versions of Sections 8 and 9: the forward code runs in the same order, the adjoint reverses it completely.
G_b 里，再把梯度那一遍反着走，最后加回壁面分量。Figure 5.1 · The data flow of the three functions. Top: residual, which computes the gradient array G from u, then r from u and G, and finally replaces the wall rows. Middle: residual_d, the same order acting on the tangent v (the state u and its gradients G frozen). Bottom: residual_b, in reverse: zero the wall components, run the flux pass backwards with the gradients' share parked in the array G_b, then run the gradient pass backwards, and finally add the wall components back.图 5.2 逐句执行这两遍循环和最后的替换。网格上节点处的紫色箭头是梯度数组 G：第一遍里它们一条边一条边地长出来，第二遍里只读不写。下方三行数组依次是状态 u*、梯度 G（每个节点两个分量）和输出 r。
Figure 5.2 executes both loops and the final replacement statement by statement. The violet arrows at the nodes of the mesh are the gradient array G: in pass 1 they grow edge by edge, in pass 2 they are only read. The three arrays below are the state u*, the gradients G (two components per node) and the output r.
播放到第二遍的某个内部面时，留意它收集了什么：除了两侧的状态 ui、uj，还有两侧的梯度 Gi、Gj，而这两个梯度又分别来自 i 和 j 的整个模板。所以一个面的通量依赖的不是两个节点，而是 i、j 两个模板里的所有节点，下一节 Jacobian 的形状就由此而来。播放到最后，壁面节点的 8 行被换成 ui−uw：在收敛解上它们恰好为 0，而循环原本留在这 8 行里的数被丢掉了。第 7 节会把这些被丢掉的数找回来：它们的相反数正是物面吸收的通量，J=−Σi∈WRi（W 是壁面节点的集合）。
When the player reaches an interior face in pass 2, note what it gathers: besides the states ui and uj on either side, the gradients Gi and Gj, and each of those came from the whole stencil of i or of j. So one face's flux depends not on two nodes but on every node in the two stencils, and that sets the shape of the Jacobian in the next section. At the very end the 8 wall rows become ui−uw: exactly 0 at the converged solution, while the numbers the loop had left in those rows are thrown away. Section 7 recovers them: with the sign reversed, they are precisely the flux the body absorbs, J=−Σi∈WRi (W the set of wall nodes).
这一节先推出残差的精确导数：切线方程与伴随方程的矩阵正是它（第 8、9 节不形成矩阵，直接作用它）；再用 Jacobi 迭代求解 Rbc(u)=0，它只要这个矩阵的对角元。
This section first derives the exact derivative of the residual, which is the matrix of the tangent and the adjoint equation (Sections 8 and 9 apply it without forming it), and then solves Rbc(u)=0 with the Jacobi iteration, which needs only its diagonal.
把一个内部面的通量对各节点值求导。对流部分与对流页面相同，只依赖两侧的 uL=ui 与 uR=uj：
Differentiate the flux of one interior face with respect to the nodal values. The convective part is as on the advection page and depends only on uL=ui and uR=uj:
扩散部分有两条路径。两点差直接依赖 ui、uj；平均梯度 Ge 则经过最小二乘，依赖 i、j 两个模板中的每一个节点（下式中 [m=j] 在 m=j 时为 1，否则为 0）：
The diffusive part has two paths. The two-point difference depends directly on ui and uj; the averaged gradient Ge goes through least squares and depends on every node of the two stencils of i and j (below, [m=j] is 1 if m=j and 0 otherwise):
所以精确 Jacobian 的第 i 行不只伸到边邻居，还伸到邻居的邻居：壁面节点 10 个非零元、中圈 12 个、外圈 11 个，全矩阵 264 个。对流页面上每行只有 i 和它的边邻居（4、6、5 个，全矩阵 120 个）。这是替换之前的 A；求解用的是替换之后的 Abc：8 个壁面行各只剩对角线上的 1，全矩阵 192 个非零元。
So row i of the exact Jacobian reaches not only the edge neighbours but also their neighbours: 10 non-zeros for a wall node, 12 for a middle-ring node, 11 for an outer-ring node, 264 in the whole matrix. On the advection page each row held only i and its edge neighbours (4, 6 and 5 entries, 120 in all). That is A before the replacement; the solve uses Abc after it, whose 8 wall rows each keep only the 1 on the diagonal: 192 non-zeros in all.
只保留相邻节点之间的耦合，就得到精确 Jacobian 的一阶部分 A1：对流导数原样保留（对流通量本身就是一阶的，面上直接用两侧的节点值，没有重构），扩散导数只保留两点部分，梯度路径整个丢掉。很多隐式求解器就拿它作迭代矩阵。本页的求解不用它：它在这里只用来看清梯度路径添了什么（图 6.1），以及第 8 节里，拿它的转置求伴随会错多少。
Keep only the couplings between adjacent nodes and what is left is the first-order part A1 of the exact Jacobian: the convective derivatives as they are (the convective flux is itself first order: it uses the two nodal values directly, with no reconstruction), and of the diffusive derivative only the two-point part, the gradient path dropped entirely. Many implicit solvers take it as their iteration matrix. The solve on this page does not use it: here it only shows what the gradient path adds (Figure 6.1) and, in Section 8, how wrong an adjoint built on its transpose would be.
A1 的非零结构与对流页面的 Jacobian 相同，全矩阵 120 个非零元；壁面行换成 êiT 之后，A1,bc 有 96 个。
A1 has the non-zero pattern of the advection page's Jacobian, 120 non-zeros in all; with the wall rows replaced by êiT, A1,bc has 96.
第一步，一个面的通量读到哪些未知量。内部对偶面的通量读两个状态 ui、uj 和两个梯度 Gi、Gj；每个梯度又读它自己的节点和它的全部边邻居。所以 he 依赖 i、j 以及两个模板里的所有节点（图 6.2 左，对偶面 e20：7 个）。远场半面没有扩散通量，只读 ui；壁面半面被跳过。
Step 1, what one face flux reads. The flux of an interior dual face reads two states, ui and uj, and two gradients, Gi and Gj; each gradient in turn reads its own node and all its edge neighbours. So he depends on i, j and every node of the two stencils (Figure 6.2, left: the dual face e20, 7 nodes). A far-field half-face carries no diffusive flux and reads ui only; the wall half-faces are skipped.
第二步，一行是它所有面的并集。第 i 行收集 ±∂h/∂um，求和取遍节点 i 的所有面，所以它的非零列就是这些面读到的未知量的并集：i 自己、它的边邻居，以及邻居的模板成员，也就是邻居的邻居。
Step 2, a row is the union over its faces. Row i collects ±∂h/∂um over all faces of node i, so its non-zero columns are the union of what those faces read: i itself, its edge neighbours, and the members of the neighbours' stencils, the neighbours of neighbours.
第三步，按节点的种类数。每一圈 8 个节点，同一圈的行非零元个数相同：壁面节点 10 个，中间一圈 12 个，远场节点 11 个。一阶部分 A1 只保留自己和边邻居：4、6、5 个。替换壁面行时，每个壁面行只剩对角元 1：
Step 3, count by kind of node. Each ring has 8 nodes, and the rows of one ring have the same number of non-zeros: 10 for a wall node, 12 on the middle ring, 11 for a far-field node. The first-order part A1 keeps only the node and its edge neighbours: 4, 6 and 5. Replacing the wall rows leaves each of them only its diagonal 1:
这就是 Abc 的来历：非壁面的 16 行原样来自 A，8 个壁面行换成 êiT（图 6.2 右）。列不受影响：壁面节点的列仍然有非壁面行里的非零元。
That is where Abc comes from: its 16 rows off the wall are the rows of A as they are, and its 8 wall rows are replaced by êiT (Figure 6.2, right). The columns are untouched: a wall node's column keeps its non-zeros in the rows off the wall.
核对：按第一、二步从每个面读到的集合拼出的 24 行，与 jacobian() 算出的精确矩阵逐行比较，每一行的非零元个数都相同（本页代码在构建时计算）。
Check: the 24 rows assembled from what each face reads (steps 1 and 2), compared row by row with the exact matrix of jacobian(), have the same number of non-zeros in every row (computed by this page's code at build time).
图 6.3 把这种结构动起来。前向：给一个节点一个切线种子（v=êj），第一遍只点亮碰到它的边，把差值分给两端的梯度，第二遍再点亮所有碰到这个状态或被改变梯度的对偶面，写进它们两侧的行，最后壁面行被 vi 覆盖，得到 Abcêj，即 Abc 的第 j 列。反向：给第 i 个残差一个伴随种子（w=êi），循环前清零壁面分量，第二遍倒过来只读第 i 行的面，第一遍倒过来再经最小二乘权重读到模板里的每个节点，最后加回壁面分量，得到 AbcTêi，即 Abc 的第 i 行。默认选 n9，它的行正是 图 6.1 左边那一行。替换壁面行使行与列连非零位置都不再相同：n9 的行有 12 个非零元，列只有 9 个，列里落在壁面行上的 3 个被替换清掉了；壁面节点 n1 的行只剩对角线上的 1，它的列却还有 6 个非零元。
Figure 6.3 sets this structure in motion. Forward: seed one node with a tangent (v=êj); pass 1 lights up only the edges that touch it and hands the difference to the gradients at both ends, pass 2 then lights up every dual face that touches the seeded state or a changed gradient and writes into the rows on both of its sides, the wall rows are finally overwritten by vi, and the result is Abcêj, column j of Abc. Reverse: give residual i an adjoint seed (w=êi); the wall components are zeroed before the loops, pass 2 reversed reads only the faces of row i, pass 1 reversed then reaches, through the least-squares weights, every node of the stencils, the wall components are added back, and the result is AbcTêi, row i of Abc. The default is n9, whose row is exactly the one on the left of Figure 6.1. With the wall rows replaced, a row and a column no longer even share their non-zero positions: row 9 has 12 non-zeros and column 9 only 9, the replacement having cleared the 3 that fell on wall rows; the row of wall node n1 keeps only the 1 on its diagonal, while its column still has 6 non-zeros.
本页的 primal、切线、伴随三个方程（后两个在第 10、11 节）用的是同一个 Jacobi 迭代：每一轮算一次残差，再逐节点除以精确 Jacobian 的对角元：
On this page the primal, the tangent and the adjoint equation (the last two in Sections 10 and 11) are all solved by one and the same Jacobi iteration: every sweep evaluates the residual once and divides it, node by node, by the diagonal of the exact Jacobian:
diag(Abc) 不需要把 Abc 造出来：diagonal(u) 就是残差的那个边循环，每个对偶面只留下落在对角线上的两项，即对流导数、两点部分，以及经由最小二乘梯度的 −12νk·(∂Gi/∂ui+∂Gj/∂ui)。壁面行是替换行，(Abc)ii=1，残差是 ui−uw，一轮就把 uw 放上去；代码直接写入 uw，使 ui=uw 逐位成立。所以一轮只有一次 residual 和一次 diagonal：不形成矩阵，也不解线性方程组。初值取 u≡1（壁面节点取 uw），判据是残差的均方根 ‖Rbc‖<10−13（图 6.4）。
diag(Abc) does not need Abc to be formed: diagonal(u) is the residual's own edge loop, keeping of each dual face only the two entries that land on the diagonal: the convective derivative, the two-point part, and the term through the least-squares gradients, −12νk·(∂Gi/∂ui+∂Gj/∂ui). A wall row is a replaced row, with (Abc)ii=1 and residual ui−uw, so one sweep puts uw there; the code writes uw in directly, so that ui=uw holds to the bit. So a sweep is one residual and one diagonal: no matrix is formed and no linear system is solved. The initial guess is u≡1 (uw on the wall nodes), and the criterion is the rms of the residual, ‖Rbc‖<10−13 (Figure 6.4).
为什么用 Jacobi，而不用 LU。本页只有 24 个未知量，理论上最快的办法是直接分解：把精确的 Abc 整个造出来，用 Newton 法，每一步做一次稠密 LU 分解，k=4 时残差就已低于 10−13。可那样就说明不了 matrix-free。Jacobi 是最简单的迭代法：对矩阵它只要对角元，其余的一切都经由残差作用进来（切线方程里是 residual_d，伴随方程里是 residual_b），这正是 matrix-free 的意思。代价是线性收敛：在解附近，误差每一轮乘以迭代矩阵 I−diag(Abc)−1Abc，速率由它的谱半径决定，本算例是 0.787，判据在 k=121 时满足。切线与伴随的迭代矩阵与它同谱，所以进入渐近阶段以后，三条收敛曲线以同一个速率下降（第 11 节）。
Why Jacobi, and not LU. With only 24 unknowns the fastest method in principle is a direct factorisation: form the whole exact Abc, run Newton, and factorise with one dense LU per step; the residual is below 10−13 by k=4. But that would show nothing about matrix-free. Jacobi is the simplest iterative method there is: of the matrix it needs only the diagonal, and everything else comes in through the residual (through residual_d in the tangent equation and residual_b in the adjoint equation), which is exactly what matrix-free means. The price is linear convergence: near the solution the error is multiplied every sweep by the iteration matrix I−diag(Abc)−1Abc, at a rate set by its spectral radius, 0.787 here, and the criterion is met at k=121. The tangent and adjoint iteration matrices have the same spectrum, so once in the asymptotic phase all three convergence curves fall at the same rate (Section 11).
对角元里的梯度路径：本页可有可无，格心页面上不可少。只取一阶部分 A1,bc 的对角元，迭代照样收敛到同一个 u*（残差是精确的），只是速率变了。在这里几乎看不出来：梯度路径最多只改变一个对角元的 3.9%，谱半径是 0.786 而不是 0.787。格心页面上它能改变一个对角元的 15.3%，谱半径从 0.915 升到 0.990，要多跑九倍左右的轮次。
The gradient path in the diagonal: optional here, essential on the cell-centred page. Take the diagonal of the first-order part A1,bc instead and the iteration still converges to the same u* (the residual is exact), only at another rate. Here that hardly shows: the gradient path changes a diagonal entry by at most 3.9%, and the spectral radius is 0.786 instead of 0.787. On the cell-centred page it changes one by up to 15.3%, and the spectral radius rises from 0.915 to 0.990, about nine times as many sweeps.
固定轮次，但会检查。代码固定跑 300 轮，而不是一满足判据就退出：余下的轮次把 u* 一直带到舍入误差平台（从 k=162 起），页面上每一个实测常数都是在这些轮次下量出来的；得到的 u* 与 Newton 法的解最多差 1.1×10−16。它记录判据在第几轮满足，若 300 轮内从未满足就直接报错，而不会把一个未收敛的状态当成 u* 返回。
A fixed sweep count, but a checked one. The code runs a fixed 300 sweeps rather than exiting at the criterion: the remaining sweeps take u* down to its round-off floor (from k=162 on), and every measured constant on this page was produced under those sweeps; the u* it returns is within 1.1×10−16 of Newton's. It records the sweep at which the criterion was met, and raises an error if 300 sweeps never meet it, rather than returning an unconverged state as u*.
| 初值处的 ‖Rbc‖ | ‖Rbc‖ at the initial guess | 0.575 |
|---|---|---|
| 谱半径 ρ(I − diag(Abc)−1Abc) | spectral radius ρ(I − diag(Abc)−1Abc) | 0.787 |
| 满足 ‖Rbc‖ < 10⁻¹³ 的轮次 | sweep meeting ‖Rbc‖ < 10⁻¹³ | k = 121（Newton 加稠密 LU：4）k = 121 (Newton with a dense LU: 4) |
| 舍入误差平台的起点 | where the round-off floor starts | k = 162 |
| 300 轮之后的 ‖Rbc‖ | ‖Rbc‖ after the 300 sweeps | 7.0e−17 |
| 壁面节点上的 u* | u* on the wall nodes | 0 （逐位）(to the bit) |
| 其余节点上的 u* | u* on the other nodes | 0.151 – 0.808 |
第一步，把一轮写成一个映射：uk+1=uk−diag(Abc(uk))−1Rbc(uk)。diag(Abc) 随 u 变，因为精确 Jacobian 随 u 变（对流导数 cL、cR 里有 u）。
Step 1, write one sweep as a map: uk+1=uk−diag(Abc(uk))−1Rbc(uk). diag(Abc) changes with u, because the exact Jacobian does (the convective derivatives cL, cR contain u).
第二步，在收敛解附近线性化。记第 k 轮的误差为 ek=uk−u*，并用 Rbc(u*)=0：
Step 2, linearise near the converged solution. Write ek=uk−u* for the error after k sweeps and use Rbc(u*)=0:
diag(Abc) 的变化乘的是 Rbc(u*)=0，所以在一次项里消失了：收敛解附近，一轮就是把误差乘以迭代矩阵 I−diag(Abc)−1Abc。
The change of diag(Abc) multiplies Rbc(u*)=0 and drops out of the linear term: near the converged solution one sweep multiplies the error by the iteration matrix I−diag(Abc)−1Abc.
第三步，重复 k 轮。把初始误差按迭代矩阵的特征向量展开，每一轮每个分量乘上它的特征值 μ；模最大的那个最后留下来：
Step 3, repeat for k sweeps. Expand the initial error in the eigenvectors of the iteration matrix; every sweep multiplies each component by its eigenvalue μ, and the one of largest modulus is what remains:
所以（局部）收敛的条件是谱半径 ρ<1，而 ρ 就是渐近段每一轮的收缩因子。
So the (local) condition for convergence is a spectral radius ρ<1, and ρ is the factor by which each sweep shrinks the error in the asymptotic phase.
第四步，一个充分条件，以及它为什么在这里不成立。迭代矩阵的对角元为零，非对角元是 −(Abc)ij/(Abc)ii。Gershgorin 圆盘定理说，每个特征值都落在某个以 0 为心、半径为 σi 的圆盘里：
Step 4, a sufficient condition, and why it fails here. The iteration matrix has a zero diagonal and off-diagonal entries −(Abc)ij/(Abc)ii. The Gershgorin disc theorem puts every eigenvalue inside one of the discs around 0 of radius σi:
若每一行都严格对角占优（σi<1），收敛就有保证。本页却有 11 行 σi>1，最大 1.113：这个界说明不了什么。收敛要看特征值本身（图 6.5）：它们全在单位圆内，ρ=0.7874，第二大的模是 0.6523。8 个壁面行的迭代矩阵行全为零（(Abc)ii=1，行是 êiT），贡献 8 个零特征值：一轮就把壁面值放到位。
If every row were strictly diagonally dominant (σi<1), convergence would be guaranteed. Here 11 rows have σi>1, the largest 1.113: the bound says nothing. Convergence rests on the eigenvalues themselves (Figure 6.5): all of them lie inside the unit circle, ρ=0.7874, and the second largest modulus is 0.6523. The 8 wall rows of the iteration matrix are zero ((Abc)ii=1 and the row is êiT), which gives 8 zero eigenvalues: one sweep puts the wall values in place.
第五步，对角元里的梯度路径。改用一阶部分的对角元 diag(A1,bc)，迭代矩阵变成 I−diag(A1,bc)−1Abc。本页梯度那一项最多只占对角元的 3.9%（图 6.6），谱半径从 0.7874 变成 0.7864，几乎不变。
Step 5, the gradient path in the diagonal. With the first-order part's diagonal diag(A1,bc) the iteration matrix becomes I−diag(A1,bc)−1Abc. On this page the gradient term is at most 3.9% of a diagonal entry (Figure 6.6), and the spectral radius goes from 0.7874 to 0.7864, hardly a change.
核对：在 u* 附近对一轮 Jacobi 映射取中心差分（步长 10−6），与 I−diag(Abc)−1Abc 最多差 5.6×10−11；由特征值算出的 ρ=0.787412，迭代在第 60 到 100 轮实测的收缩因子是 0.787412（本页代码在构建时计算）。
Check: central differences of one Jacobi sweep near u* (step 10−6) differ from I−diag(Abc)−1Abc by at most 5.6×10−11; the eigenvalues give ρ=0.787412, and the iteration shrinks by 0.787412 per sweep measured over k=60 to 100 (computed by this page's code at build time).
物面吸收的总通量是这类问题最常见的目标函数：传热问题里它是热流，流动问题里对应表面摩擦。强施加恰好把它藏了起来，因为壁面节点的守恒方程已经换掉了。但在换掉之前，边循环已经把这一行的通量和算好了（Eint 是内部对偶面的集合）：
The total flux the body absorbs is the most common objective of such problems: the heat flux in heat transfer, the counterpart of skin friction in a flow. Strong imposition hides it, because the wall nodes' conservation equations have been swapped out. But before the swap, the edge loop had already computed each of those rows' flux sums (Eint is the set of interior dual faces):
而这个控制体真正的平衡还要加上两个壁面半面上的通量 qi（从 Vi 流进物面）：Ri(u)+qi=0。所以
while the true balance of that control volume also includes the flux qi through its two wall half-faces (out of Vi, into the body): Ri(u)+qi=0. Hence
这称为反力形式：物面要吸收多少，才能让每个壁面节点自己的守恒方程成立。它与离散格式一致：把 24 个控制体的真实平衡加起来，内部对偶面的通量两两抵消，剩下的只有 J 与远场半面通量之和，所以 J 恰好等于从远场流入的净通量。实测 J=2.748718，远场半面的通量之和是 −2.748718，两者相加为 0。
This is the reaction form: how much the body must absorb for every wall node's own conservation equation to hold. It is consistent with the discrete scheme: add up the true balances of all 24 control volumes, the fluxes through interior dual faces cancel in pairs, and all that remains is J plus the sum over the far-field half-faces, so J is exactly the net flux entering through the far field. Measured: J=2.748718, the far-field half-faces sum to −2.748718, and the two add up to 0.
另一种常见做法是用单侧差分直接估计壁面法向导数。每个壁面节点 i 只有一个不在壁面上的边邻居 i′（中圈上的那个节点），于是
A common alternative estimates the wall-normal derivative directly with a one-sided difference. Each wall node i has exactly one edge neighbour i′ off the wall (its node on the middle ring), and
其中 ℓi 是节点 i 两个壁面半面的长度之和。在本网格上 J1=0.919528，只有一致值的 33%。原因有二：壁面第一层网格很粗，单侧差分的截断误差大；更根本的是，它只看壁面节点的一个邻居，而离散格式里流进物面的通量是由整个壁面控制体的平衡决定的，这个控制体从其他面上流进来的东西（包括被对流带进来的）它都没有计入。图 7.1 把两者逐个壁面节点并排画出。
where ℓi is the combined length of node i's two wall half-faces. On this mesh J1=0.919528, only 33% of the consistent value. There are two reasons: the first layer of cells at the wall is coarse, so the one-sided difference has a large truncation error; and, more fundamentally, it looks at one neighbour of the wall node, whereas in the discrete scheme the flux into the body is set by the balance of the whole wall control volume, and whatever enters that volume through its other faces (convection included) is left out. Figure 7.1 shows the two side by side, wall node by wall node.
把网格一致加密（第 15 节的网格族），两者从两侧互相靠拢：反力通量依次为 2.749、2.454、2.291、2.201、2.153、2.127，单侧公式为 0.920、1.314、1.625、1.821、1.925、1.974（24、80、288、1088、4224、16640 个节点）。反力通量在每一层都与远场通量精确平衡；单侧公式在任何一层都不满足这一点。
Refine the mesh uniformly (the family of Section 15) and the two approach each other from opposite sides: the reaction flux runs 2.749, 2.454, 2.291, 2.201, 2.153 and 2.127, the one-sided formula 0.920, 1.314, 1.625, 1.821, 1.925 and 1.974 (24, 80, 288, 1088, 4224 and 16640 nodes). The reaction flux balances the far-field flux exactly on every level; the one-sided formula does so on none.
J=−Σi∈WRi(u) 是残差分量的线性组合，所以 gT=∂J/∂u 就是替换之前的 A 的 8 个壁面行之和，再取负：
J=−Σi∈WRi(u) is a linear combination of residual components, so gT=∂J/∂u is minus the sum of the 8 wall rows of A before the replacement:
在代码里它就是一次伴随循环：壁面节点的权重取 −1、其余取 0，走替换之前的 residual_b。精确的 A 每行伸到邻居的邻居，所以 g 在全部 24 个节点上都非零，而单侧公式的导数只在 8 个中圈节点 i′ 上非零。
In code it is one adjoint loop: weights −1 on the wall nodes and 0 elsewhere, through residual_b before any replacement. Each row of the exact A reaches the neighbours of neighbours, so g is non-zero on all 24 nodes, whereas the one-sided formula's derivative is non-zero only on the 8 middle-ring nodes i′.
第一步，目标函数是替换之前壁面行的线性组合：J=−Σi∈WRi(u,X)。对状态求偏导（X 不动），每一项就是 A 的一行：
Step 1, the objective is a linear combination of the wall rows before the replacement: J=−Σi∈WRi(u,X). With X held fixed, each term differentiates into one row of A:
第二步，哪些分量非零。A 的每一行都伸到邻居的邻居（上一节），8 个壁面行合起来覆盖了全部 24 个节点：g 处处非零（图 7.2 左）。单侧公式 J1=Σi∈Wν(ui′−uw)ℓi/|xi′−xi| 只读 8 个中间节点，它的导数 ∂J1/∂ui′=νℓi/|xi′−xi| 也只在那里非零（右）。
Step 2, which components are non-zero. Every row of A reaches the neighbours of neighbours (Section 6), and the 8 wall rows together cover all 24 nodes: g is non-zero everywhere (Figure 7.2, left). The one-sided formula J1=Σi∈Wν(ui′−uw)ℓi/|xi′−xi| reads only the 8 middle-ring nodes, and its derivative ∂J1/∂ui′=νℓi/|xi′−xi| is non-zero only there (right).
第三步，代码怎么算。gT=wTA 取 w=−1W，正是替换之前的 residual_b（第 9 节，不清零也不加回）作用在这个 w 上：每个面通量的伴随 h̄=wi−wj 只在壁面节点与它的中圈邻居之间的 8 个对偶面上是 −1，其余的面上都是 0；两个壁面节点之间的 8 个面上两个 −1 相消，因为这种通量从一个壁面行流出、流进另一个壁面行。
Step 3, how the code computes it. gT=wTA with w=−1W is exactly residual_b before any replacement (Section 9, no zeroing, no adding back) applied to that w: each face flux's adjoint h̄=wi−wj is −1 on the 8 dual faces between a wall node and its middle-ring neighbour and 0 on every other face: on the 8 faces between two wall nodes the two −1 cancel, because such a flux leaves one wall row and enters the other.
第四步，对坐标的偏导（状态不动）。同一个 J=−Σi∈WRi 也直接依赖坐标：
Step 4, the partial derivative with respect to the coordinates, the state held fixed. The same J=−Σi∈WRi depends on the coordinates directly as well:
它是一次对坐标的反向扫描，种子就是同一个 w=−1W 给出的 h̄：第 12 节的 φ 在壁面上取 −1，正是这一项。
It is one reverse sweep over the coordinates, seeded by the h̄ of the same w=−1W: that is exactly the −1 the φ of Section 12 takes on the wall.
核对：按上面的公式逐项算出的 g 与代码 dJdu()（即 residual_b）最多差 4.4×10−16，与对 J 的中心差分（步长 10−6）最多差 2.3×10−10；∂J/∂X 的反向扫描与状态不动时移动坐标的中心差分最多差 4.2×10−10，它的分量最大为 0.624（本页代码在构建时计算）。
Check: the g computed term by term from the formulas above differs from the code's dJdu() (that is, residual_b) by at most 4.4×10−16, and from central differences of J (step 10−6) by at most 2.3×10−10; the reverse sweep for ∂J/∂X differs from central differences that move the coordinates with the state held fixed by at most 4.2×10−10, whose largest component is 0.624 (computed by this page's code at build time).
不能只对连续积分写一个形式导数就结束。先确定代码实际使用的求积点、权重、重构状态及边界值，再对这个标量表达式求导。一般的求积和及其链式法则是：
A formal derivative of a continuous integral is not enough. First identify the quadrature points, weights, reconstructed states and boundary values actually used by the code, then differentiate that scalar expression. A general quadrature sum and its chain rule are:
本页的反力形式 是替换前残差行的线性组合：，（替换之前的第 行）， 取遍壁面节点。
The reaction form of this page is a linear combination of residual rows before replacement: and (row before replacement), over the wall nodes.
求 时固定所有设计变量；求 时固定全部离散状态槽位。固定解析目标函数并不一定意味着移动网格上的目标采样值不变：如果采样点随设计移动，还要保留该位置依赖；固定网格及固定采样值时，相应项才为零。目标中的体积权重已经进入 ，不应在伴随右端再乘一次体积。
Hold every design variable fixed for , and every stored state entry fixed for . A fixed analytic target function need not have fixed samples on a moving mesh: if the evaluation locations move with the design, keep that dependence; the corresponding terms vanish only when the mesh and samples are fixed. Objective volume weights already enter and must not be multiplied into the adjoint right-hand side again.
伴随右端本身也可以由一次反向扫描得到：以 为种子，把上面的求积和逐项反过来执行，
The adjoint right-hand side is itself obtained by a reverse sweep: seed and run the terms of the quadrature sum above backwards,
从 、 开始，扫完全部求积点后 、：一个状态被多个求积点读取，所以要先置零、再累加。前向的 与反向的这两个结果满足同一个对偶关系 。
Starting from and , after all quadrature points and : a state is read by several quadrature points, so its adjoint is zeroed first and then accumulated. The forward and these two reverse results satisfy one duality, .
以设计变量 D 中的网格节点坐标 X 为例，共 48 个分量，是第 1 节那 81 个分量的一部分；其余分量的推导相同，第 10 节末尾逐类给出。状态由 Rbc(u,X)=0 决定，目标函数 J(u,X) 对 X 的全导数是
Take the mesh node coordinates X of the design variables D as the example: 48 components, part of the 81 of Section 1; the other components go the same way, kind by kind at the end of Section 10. The state is fixed by Rbc(u,X)=0, and the total derivative of the objective J(u,X) with respect to X is
推导与对流页面第 8 节相同。这里要强调的是：伴随方程里的矩阵是精确的 Abc，不是它的一阶部分 A1,bc（第 6 节；很多隐式求解器留用的就是它）。拿一阶部分的转置 A1,bcT 解出来的 ψ 满足的是另一个方程：本算例里它与精确的 ψ 最多差 0.033，给出的形状梯度最多差 2.6e−3（最大分量的 0.8%），而精确伴随与有限差分的差在 10−9 量级（第 13 节）。精确矩阵每行伸到邻居的邻居，形成它既麻烦又浪费，所以和对流页面一样，只需要能算乘积 Abcv 与 AbcTw 的两个函数：residual_d 与 residual_b。这一节是前者。
The derivation is that of Section 8 of the advection page. The point to stress here is that the matrix in the adjoint equation is the exact Abc, not its first-order part A1,bc (Section 6), which many implicit solvers keep. A ψ solved with the transpose of the first-order part, A1,bcT, satisfies a different equation: here it misses the exact ψ by up to 0.033, and its shape gradient misses by up to 2.6e−3 (0.8% of the largest component), whereas the exact adjoint agrees with finite differences to about 10−9 (Section 13). The exact matrix reaches the neighbours of neighbours in every row, awkward and wasteful to form, so, as on the advection page, all that is needed are two functions that compute the products Abcv and AbcTw: residual_d and residual_b. This section is the first of them.
记号。本页的前向与伴随都是对离散代码求导，用自动微分的记号：程序变量 x 的切线（前向模式的导数）记作 ẋ，伴随（反向模式的导数）记作 x̄；伴随先清零（x̄←0），再按语句的逆序累加（x̄+=…）。代码与演示里用 Tapenade 的名字 xd≡ẋ、xb≡x̄：residual_d 的输入 v 就是状态的切线 u̇，输出 r_d 就是 ṙ=Abcu̇。δ 只表示连续函数的变分，只有附录 A 的连续伴随用变分符号。
Notation. The forward and the adjoint mode of this page differentiate the discrete code, and they are written in the notation of algorithmic differentiation: the tangent (forward-mode derivative) of a program variable x is ẋ, its adjoint (reverse-mode derivative) x̄; an adjoint is cleared first (x̄←0) and then accumulated (x̄+=…) in the reverse order of the statements. The code and the players use Tapenade's names, xd≡ẋ and xb≡x̄: the input v of residual_d is the state's tangent u̇, and its output r_d is ṙ=Abcu̇. The symbol δ is kept for the variation of a continuous function: only the continuous adjoint of Appendix A uses the variation symbol.
沿残差的执行顺序，逐遍线性化。第一遍的梯度对 u 是线性的，所以它的切线就是把同一个循环、同一组权重作用在 u̇ 上。第二遍的通量对两侧状态和平均梯度 Ge 求导，对流部分的系数 cL、cR 与对流页面相同。最后，替换行 ui−uw 的切线就是 u̇i：
Linearise the residual pass by pass, in the order it runs. The gradients of pass 1 are linear in u, so their tangent is the same loop with the same weights, applied to u̇. The flux of pass 2 is differentiated with respect to the two states and the averaged gradient Ge, with the convective coefficients cL, cR of the advection page. Finally, the tangent of a replaced row ui−uw is simply u̇i:
每条从 u 到 r 的路径都被求了导，包括经过梯度的那一条，所以这就是精确的 ṙ=Abcu̇，即 Abcv。去掉 Ġe·k 一项，得到的就是一阶部分 A1,bcv。两者的差正是一阶部分丢掉的东西：取 v=ê9，差的最大值为 0.115（Abcê9 本身最大 4.03），非零分量落在 n9 两圈之内的 9 个节点上。
Every path from u to r is differentiated, the one through the gradients included, so this is the exact ṙ=Abcu̇, that is Abcv. Drop the Ġe·k term and what is left is the first-order part A1,bcv. The difference between the two is exactly what the first-order part drops: with v=ê9 it reaches 0.115 (against 4.03 for Abcê9 itself), on 9 nodes within two rings of n9.
图 8.1 逐句执行 residual_d。三行数组依次是切线 v（即 u̇）、它的梯度 G_d（即 Ġ）和输出 r_d（即 ṙ）。点击控制体可以改种子；选壁面节点的 ê1 时，第二遍照样把 v1 带进邻居的行，最后第 1 行被覆盖成 v1：替换行只改自己那一行，不改那一列。矩阵视图里的显式 Abc 只为校对而形成。
Figure 8.1 executes residual_d statement by statement. The three arrays are the tangent v (u̇), its gradients G_d (Ġ) and the output r_d (ṙ). Click a control volume to change the seed; with the wall node's ê1, pass 2 still carries v1 into its neighbours' rows, and at the end row 1 is overwritten by v1: replacing a row changes that row only, not the column. The explicit Abc in the matrix view is formed only for the check.
由链式法则，residual_d(u,v) 就是 ddtRbc(u+tv) 在 t=0 处的值，也就是 Abcv。它与显式乘积的差在舍入量级（播放到最后，矩阵视图下方给出两列的最大差），而全程没有形成矩阵。
By the chain rule, residual_d(u,v) is ddtRbc(u+tv) at t=0, that is Abcv. It differs from the explicit product only by round-off (play to the end: the matrix view reports the largest difference between the two columns), and no matrix was ever formed.
上面的切线公式已经够算 Abcu̇。第 10 节的右端项和第 12 节的设计导数还要更多：每一种边界面，以及设计输入（来流值、源项、壁值、坐标）进入残差的每一处。下面把 residual 的每一句（第 5 节的伪代码）单独求导。只用一条规则，就是链式法则：一句的每个输入 x 带着自己的切线 ẋ，输出的切线是各输入的切线乘这一句的偏导数再相加。
The tangent formulas above are enough for Abcu̇. The right-hand sides of Section 10 and the design derivatives of Section 12 need more: every kind of boundary face, and every place where a design input (the free-stream value, a source, a wall value, a coordinate) enters the residual. Below, each statement of residual (the pseudocode of Section 5) is differentiated on its own. One rule does it all, the chain rule: every input x of a statement carries its tangent ẋ, and the output's tangent is the sum of the input tangents times the statement's partial derivatives.
第一步，第一遍（第 4 至 6 行）：沿边 (i,j) 的差值 Δ=uj−ui，再累加 Gi+=wijΔ 与 Gj+=−wjiΔ。输入是两个状态和两个权重，所以
Step 1, pass 1 (lines 4 to 6): the difference along edge (i,j), Δ=uj−ui, and the two accumulations Gi+=wijΔ and Gj+=−wjiΔ. The inputs are two states and two weights, so
权重只依赖几何，ẇ 只在坐标带切线时不为零（第 12 节推出它）。只对状态求切线时只剩各式的第一项：第一遍的切线就是同一个循环作用在 u̇ 上，即 residual_d 的第 4 至 6 行。
The weights depend on the geometry only, so ẇ is non-zero only for a coordinate seed (Section 12 derives it). For a state tangent only the first terms remain: pass 1's tangent is the same loop acting on u̇, lines 4 to 6 of residual_d.
第二步，第二遍的收集（第 10、11 行）：取两侧状态与平均梯度，切线照样取。右状态随面的种类而定：内部对偶面上 u̇R=u̇j，远场入流半面上 u̇R=u̇∞，出流半面上 uR=uL，所以 u̇R=u̇L。平均梯度 Ġe=12(Ġi+Ġj)。壁面半面被跳过。对状态求切线时 u∞ 是常数，u̇∞=0。
Step 2, the gathers of pass 2 (lines 10 and 11): the two states and the averaged gradient, and their tangents likewise. The right state depends on the kind of face: on an interior dual face u̇R=u̇j, on a far-field inflow half-face u̇R=u̇∞, and on an outflow half-face uR=uL, so there u̇R=u̇L. The averaged gradient gives Ġe=12(Ġi+Ġj). The wall half-faces are skipped. In the state's tangent u∞ is a constant, u̇∞=0.
第三步，对流通量（第 12 行的前一半）。vn=β·n 与 |n| 都随 n 变：v̇n=β·ṅ，ddt|n|=n·ṅ/|n|。于是
Step 3, the convective flux (the first half of line 12). Both vn=β·n and |n| move with n: v̇n=β·ṅ and ddt|n|=n·ṅ/|n|. Hence
远场半面只是这条式子的特例（壁面半面被跳过）。出流半面上 uR=uL（第 4 节），两个槽位合成 (cL+cR)u̇L=vnuLu̇L：代码在出流半面上的 cL 就是这个合并后的系数。入流半面上的 cRu̇∞ 是 u∞ 进入残差的唯一一处，它就是第 10 节的 f。
The far-field half-faces are special cases of this formula (the wall half-faces are skipped). On an outflow half-face uR=uL (Section 4), and the two slots merge into (cL+cR)u̇L=vnuLu̇L: the code's cL on an outflow half-face is this merged coefficient. On an inflow half-face the term cRu̇∞ is the one place where u∞ enters the residual: it is the f of Section 10.
第四步，扩散通量（第 12 行的后一半）。hv 读四样东西：两侧状态、平均梯度、p 和 k；而 p=ξ·n/|ξ|2、k=n−pξ 由 n 与 ξ 算出（第 3 节）：
Step 4, the diffusive flux (the second half of line 12). hv reads four things: the two states, the averaged gradient, p and k; and p=ξ·n/|ξ|2 and k=n−pξ come from n and ξ (Section 3):
把后两式代进前一式，按 ṅ 与 ξ̇ 整理，几何部分正是第 12 节的 ∂hv∂n·ṅ+∂hv∂ξ·ξ̇。远场半面上 hv=0，切线也是 0。
Put the last two into the first and collect by ṅ and ξ̇: the geometric part is exactly Section 12's ∂hv∂n·ṅ+∂hv∂ξ·ξ̇. On a far-field half-face hv=0, and so is its tangent.
第五步，分发（第 13、14 行）、源项与壁面行。累加是线性的，切线照样累加；源项 si（第 11 节的设计变量）进入 Ri←Ri−si；壁面行（第 16 行）ri←ui−uw,i 最后执行：
Step 5, the scatters (lines 13 and 14), the source and the wall rows. Accumulations are linear, and their tangents accumulate the same way; a source si (a design variable of Section 11) enters as Ri←Ri−si; and the wall rows (line 16), ri←ui−uw,i, run last:
第六步，哪个种子喂哪一句。把第一至五步按输入读一遍，就知道每一类设计分量的右端项 f=∂Rbc/∂Dk 从哪几句来（第 10 节末尾逐个解它们）：
Step 6, which seed feeds which statement. Read steps 1 to 5 by their inputs and it is clear which statements the right-hand side f=∂Rbc/∂Dk of each kind of design component comes from (the end of Section 10 solves them one by one):
图 8.2 把 residual_d 整个跑一遍，种子取第 10 节切线方程的解 v=du/du∞：第一遍逐边点亮梯度的切线，第二遍逐面点亮通量的切线，右边的柱是输出 ṙ 在每个节点上的累计值。v 满足 Abcv=−f，所以跑完以后只有 5 个入流节点上 ṙ=−f 不为零，其余的正负抵消：这一次前向扫描也就核对了第 10 节的解。
Figure 8.2 runs residual_d once through, seeded with the solution v=du/du∞ of Section 10's tangent equation: pass 1 lights up the gradients' tangents edge by edge, pass 2 the fluxes' tangents face by face, and the bars on the right are the output ṙ accumulated in each node. Since v satisfies Abcv=−f, at the end ṙ=−f is non-zero only on the 5 inflow nodes and everything else has cancelled: one forward sweep that also checks Section 10's solution.
核对：第一至五步逐句写成程序，与本页代码比较。状态种子下与 residual_d 的相对差最多 1.8×10−16；来流值种子下与第 10 节的 f 逐位相同；坐标种子下与第 14 节的复数步长（h=10−30）相对差最多 4.0×10−16，最小二乘权重的切线最多差 1.3×10−15；源项种子与中心差分逐位相同；动画最后一帧与 −f 最多差 4.4×10−16（本页代码在构建时计算）。
Check: steps 1 to 5, written as a program statement by statement, against this page's code. With a state seed they differ from residual_d by at most 1.8×10−16 relative; with the free-stream seed they and Section 10's f agree to the last bit; with a coordinate seed they differ from the complex step of Section 14 (h=10−30) by at most 4.0×10−16 relative, and the least-squares weights' tangents by at most 1.3×10−15; with a source seed they and central differences agree to the last bit; and the animation's last frame differs from −f by at most 4.4×10−16 (computed by this page's code at build time).
第一步，对角元只管速度，不管答案。伴随的 Jacobi 迭代 ψ←ψ−diag(Abc)−1(AbcTψ+g) 不论对角阵取什么（只要可逆），停下的地方都满足：
Step 1, the diagonal decides the speed, not the answer. Whatever the diagonal is (as long as it is invertible), the adjoint Jacobi iteration ψ←ψ−diag(Abc)−1(AbcTψ+g) can only stop where
答案由残差里的算子决定，与对角阵无关。用精确的 residual_b，但把 diag(Abc) 换成一阶部分的对角元 diag(A1,bc)，迭代 400 轮后与精确的 ψ 只差 2.2×10−16：第 6 节的 primal 也是这样，换了对角元，照样收敛到同一个 u*。
The answer is set by the operator in the residual, not by the diagonal. With the exact residual_b but the first-order part's diagonal diag(A1,bc) in place of diag(Abc), 400 sweeps end within 2.2×10−16 of the exact ψ: just as the primal of Section 6 reaches the same u* with another diagonal.
第二步，换掉算子，不动点就跟着换。若残差里用一阶部分，ψ̃ 解的是 A1,bcTψ̃=−g。减去 AbcTψ=−g，再在左边加减 A1,bcTψ：
Step 2, change the operator and the fixed point changes with it. With the first-order part in the residual, ψ̃ solves A1,bcTψ̃=−g. Subtract AbcTψ=−g and add and subtract A1,bcTψ on the left:
驱动误差的正是一阶部分丢掉的那一块，即梯度路径 Abc−A1,bc；迭代多少轮、判据多严都去不掉它。
The error is driven by exactly what the first-order part drops, the gradient path Abc−A1,bc; no number of sweeps and no tolerance removes it.
第三步，形状梯度继承这个误差。第 12 节的梯度是一次反向扫描，对种子是线性的，所以用 ψ̃ 算出的梯度与真值之差是：
Step 3, the shape gradient inherits the error. The gradient of Section 12 is one reverse sweep, linear in its seeds, so the gradient computed with ψ̃ misses the true one by
（φ̃ 与 φ 在壁面上都取 −1，所以只有非壁面节点上的差起作用。）
(φ̃ and φ both take −1 on the wall, so only the difference off the wall acts.)
核对：本页 ψ̃ 与 ψ 最多差 0.033（在壁面节点 n2 上，它不进入 dJ/dX；图里画的壁面以外的节点上最多差 0.0092），第二步的等式两边最多差 1.2×10−15；用 ψ̃ 算出的形状梯度最多差 2.6×10−3，是最大分量的 0.8%（图 8.3，本页代码在构建时计算）。
Check: on this page ψ̃ and ψ differ by at most 0.033 (at the wall node n2, which does not enter dJ/dX; on the nodes off the wall, the ones the figure draws, by at most 0.0092), and the two sides of the equation of step 2 by at most 1.2×10−15; the shape gradient computed with ψ̃ misses by up to 2.6×10−3, 0.8% of the largest component (Figure 8.3; computed by this page's code at build time).
固定网格和当前可微分支，把一个面的所有输入记为 ；它包括重构状态以及该通量需要的梯度、几何和边界数据。令 、 分别是 对存储状态和设计变量的导数， 是局部通量 Jacobian。于是每个残差行的导数按同一个有向面关联矩阵 装配，这里的求和包括相应的边界通量映射（本页的这些求和是面循环的，即覆盖之前的残差 的）；被替换的边界行另按第 9 节处理。前向模式按计算顺序传播切向量，状态种子是 ，设计种子是 ：
At fixed mesh and differentiable branch, let collect all face inputs: reconstructed states and the gradients, geometry and boundary data required by the flux. Let and be the derivatives of with respect to the stored states and the design variables, and the local flux Jacobian. Assemble each residual derivative with the same oriented face incidence , including the boundary flux maps (on this page these sums are those of the face loop, that is of the residual before any row is overwritten); replaced boundary rows are handled in Section 9. Forward mode carries the tangents in the order of the computation, with the state seed and the design seed :
无体源项时 。面长度和法向量随节点坐标变化，它们的设计依赖属于 ，不能只对状态求导。最小二乘梯度使 连接邻居的邻居，因此精确 Jacobian 的模板比一阶对流宽；matrix-free 仍只需依次执行这些局部乘积，不必存储全局矩阵。
Set when there is no volume source. Face lengths and normals move with the node coordinates; that design dependence belongs in , so differentiating the states alone is insufficient. Least-squares gradients connect to the neighbours of neighbours, so the exact Jacobian stencil is wider than first-order convection. Matrix-free evaluation performs these local products without storing the global matrix.
上面的 、 与 来自同一段短程序：残差代码对每个面依次执行四条语句。设 是面 读取的存储位置（两侧的 、，以及重构或梯度需要时它们的邻居），把这些位置上的状态排成 ：
The , and above come from one short program: the residual code runs four statements on every face. Let be the stored locations that face reads (its two sides and , and their neighbours when a reconstruction or a gradient needs them), and stack the states stored there as :
（1）取值： 的每一行只有一个 1，从全局数组里挑出 上的状态；（2）重构：由这些状态和设计变量算出面的全部输入 ，在本页是两侧状态 、 与两侧的最小二乘梯度 、，加上面积向量 、面长 与从左点到右点的向量 （节点页是两个节点；扩散通量的两点部分与偏斜部分由它们算出；这三个几何量都是节点坐标 的函数，一阶时对状态 的这几行为零）；（3）通量：（本页把面通量 记作 ）；（4）分发： 是有向关联矩阵 的第 列，在第 行为 、第 行为 ，即 、（边界面只写一行）。这四条语句是残差代码的骨架；本页代码里的每一条语句（含边界面与源项）都在本节里逐句写出，并附各自的切向量语句。
(1) Gather: each row of holds a single 1 and picks one state of out of the global array. (2) Reconstruct: the states and the design give all inputs of the face, on this page the two states , and the least-squares gradients , of the two sides, plus the area vector , the face length and the vector from the left point to the right point (the two nodes on the node page; the two-point part and the skew part of the diffusive flux are computed from them; all three geometry quantities are functions of the node coordinates , and their rows of on the states are zero). (3) Flux: (this page writes the face flux as ). (4) Scatter: is column of the oriented incidence , in row and in row , that is and (a boundary face writes one row). These four statements are the skeleton of the residual code; every statement of this page's code, boundary faces and sources included, is written out with its tangent statement in this section.
前向模式对每条语句用一次链式法则：每条语句输出的切向量是它输入的切向量的线性组合，每个中间变量都有自己的切向量：
Forward mode applies the chain rule once per statement: the tangent of each output is a linear combination of the tangents of the statement's inputs, and every intermediate variable gets a tangent of its own:
取值与分发都是线性的，它们的切向量就是同样的取值与分发作用在切向量上；只有重构与通量两条语句带来依赖当前状态的局部导数 、 与 。把（1）代入（2），（2）代入（3），（3）代入（4），对所有面求和，再减去体源项的切向量 ，就逐条语句得到 ：
Gather and scatter are linear: their tangents are the same gather and scatter applied to the tangents. Only the reconstruct and flux statements bring local derivatives that depend on the current state: , and . Substituting (1) into (2), (2) into (3) and (3) into (4), summing over all faces and subtracting the tangent of the volume source assembles statement by statement:
上面的 就是 。所以一个面对 的贡献是四个初等 Jacobian 的乘积
The above is . One face therefore contributes to the product of four elementary Jacobians,
它从右往左作用在 上，与程序执行的顺序相同：先取值，再重构、通量，最后分发。 只有几个位置，所以 只有一小块非零，matrix-free 乘积就是按面执行这几次局部运算，不存储任何全局矩阵。图 8.4 按这四条语句逐条播放本页一个面的数值。
It acts on from right to left, in the order in which the program runs: gather first, then reconstruct and flux, scatter last. holds only a few locations, so has a small nonzero block, and the matrix-free product runs these local operations face by face without storing any global matrix. Figure 8.4 plays the four statements with the numbers of one face of this page.
取面循环里第一个两侧都不在壁面上的内部面 （节点 9 在左、10 在右）；扩散通量要用两侧的最小二乘梯度，梯度又读各自的邻居，。设计种子取 ， 是来流值 （本页前向求解的方向）；状态种子取本页前向求解对同一方向给出的解 ，于是下面每个切向量都是对 的导数。表中的数由本页代码在构建时算出：基准状态、前向求解与局部导数都用本页自己的程序。
Take the first interior face of the face loop whose two sides are off the wall, (node 9 on the left, 10 on the right); the diffusive flux uses the least-squares gradients of both sides, which read their own neighbours, . Seed the design with , where is the free-stream value (the direction of this page's tangent solve), and the states with this page's tangent solution for the same direction, ; every tangent below is then a derivative with respect to . The numbers are computed by this page's code at build time: its baseline state, its tangent solve and its local derivatives.
| 语句Statement | 量Quantity | 原始值Primal value | 切向量Tangent |
|---|---|---|---|
| （1）取值(1) gather | 0 | 0 | |
| 0 | 0 | ||
| 0.1516 | 0.1853 | ||
| 0.1509 | 0.1989 | ||
| 0.2705 | 0.3695 | ||
| 0.2772 | 0.3361 | ||
| 0.2151 | 0.2617 | ||
| 0.1802 | 0.2343 | ||
| 0.2976 | 0.4016 | ||
| 设计种子design seed | 1 | 1 | |
| （2）重构(2) reconstruct | 0.1516 | 0.1853 | |
| 0.1509 | 0.1989 | ||
| 0.07305 −0.021 | 0.08339 −0.01526 | ||
| −0.03361 0.09343 | −0.05402 0.1287 | ||
| −0.4099 0.2802 | 0 0 | ||
| 0.4965 | 0 | ||
| −0.3557 0.8865 | 0 0 | ||
| （3）通量(3) flux | 0.001733 | −0.01137 | |
| （4）分发(4) scatter | 0.001733 | −0.01137 | |
| −0.001733 | 0.01137 |
不直接进入这个面：， 全部来自状态的切向量。 只作为远场入流面的右状态进入通量，所以 只在 5 个节点上非零；几何量也是 的一部分，但这个种子不移动任何节点，所以它们的切向量都是零。
does not enter this face directly: , and comes entirely from the state tangents. enters only as the right state of the far-field inflow faces, so is nonzero in 5 nodes only; the geometry is part of too, but this seed moves no node, so its tangents are zero.
所有面与源项的这些贡献加起来是替换之前的残差的切向量 。它在保留的行上于收敛解处为零；在被替换的行上它不为零（那里模最大的一项是 ），因为残差代码在面循环之后把这些行整行覆盖，覆盖之后的切向量是约束的切向量。下面前向求解的方程是整个方程组的：
Summed over all faces and the source terms, these contributions give the tangent of the raw residual, . In the retained rows it vanishes at the converged solution; in the replaced rows it does not (its entry of largest magnitude there is ), because the residual code overwrites those rows after the face loop and the tangent after the overwrite is the constraint's. The tangent equation solved further on is the system's:
按程序顺序执行四条语句，箭头上是各自的初等 Jacobian；下方是它们的乘积与形状（行 × 列，本页这个面的尺寸）。
The four statements in program order, each arrow labelled with its elementary Jacobian; below, their product with the shapes (rows × columns, this face of this page).
种子：（ 是来流值 ），。在 上： 0， 0， 0.1853， 0.1989， 0.3695， 0.3361， 0.2617， 0.2343， 0.4016。
Seeds: ( is the free-stream value ) and . On : 0, 0, 0.1853, 0.1989, 0.3695, 0.3361, 0.2617, 0.2343, 0.4016.
（1）取值：，只挑出 上的 9 个状态的切向量（上一步列出的数）。
(1) Gather: picks the tangents of the 9 states of (the numbers of the previous step).
（2）重构：。 0.1853， 0.1989， (0.08339, −0.01526)， (−0.05402, 0.1287)。 (0, 0)， 0， (0, 0)（这个种子不移动任何节点）。： 不进入这个面。
(2) Reconstruct: . 0.1853, 0.1989, (0.08339, −0.01526), (−0.05402, 0.1287). (0, 0), 0, (0, 0) (this seed moves no node). : does not enter this face.
（3）通量： −0.01137。
(3) Flux: −0.01137.
（4）分发： −0.01137，，即加上 0.01137。
(4) Scatter: −0.01137 and , that is 0.01137 added.
严格推导伴随最稳妥的办法，是先把前向写成一串线性映射的乘积，再逐个转置。第一遍是一个 2N×N 的线性算子 Γ（N=24 是未知量个数）：Ġ=Γu̇，它的第 i 个块行是 Σswis(ês−êi)T。第二遍把状态和梯度一起映到输出：ṙ=Buu̇+BgĠ。Bu 收集每个面上直接作用在两侧状态上的系数：第 i 行在第 i、j 列分别加 cL+νp、cR−νp，第 j 行加它们的相反数。Bg 收集作用在梯度上的系数：每个面在第 i 行写 −12νkT、在第 j 行写 +12νkT，两行都写进 Gi 和 Gj 两个块。最后是替换：P 把壁面行置零，E 只在壁面行的对角线上为 1，P+E=I。于是
The safest way to derive the adjoint rigorously is to write the forward mode as a product of linear maps and transpose them one by one. Pass 1 is a 2N×N linear operator Γ (N=24 unknowns): Ġ=Γu̇, whose i-th block row is Σswis(ês−êi)T. Pass 2 maps the states and the gradients together to the output: ṙ=Buu̇+BgĠ. Bu collects the coefficients acting on the states directly: on every face, row i gets cL+νp in column i and cR−νp in column j, and row j gets their negatives. Bg collects those acting on the gradients: each face writes −12νkT into row i and +12νkT into row j, both into the two blocks Gi and Gj. Last comes the replacement: P zeroes the wall rows and E is 1 on the wall rows' diagonal only, so P+E=I. Hence
用反向模式的记号：输入 w 就是 residual_d 所算残差的伴随 r̄，输出是状态的伴随 ū=AbcTr̄（代码里的 w 与 y），ū 与 Ḡ 先清零。从右往左读，每个因子都说明了代码该做什么、按什么顺序做：
In reverse-mode notation the input w is the adjoint r̄ of the residual that residual_d differentiates, and the output is the state's adjoint ū=AbcTr̄ (the code's w and y), with ū and Ḡ cleared first. Read from right to left, each factor says what the code must do, and in which order:
w̃）：替换丢掉了循环留在这些行里的值，所以不会有伴随从这些行流回去。w̃): the replacement threw away what the loops had left in those rows, so no adjoint flows back through them.图 9.1 把三个函数逐项并排。顺序是整件事的关键。A=Bu+BgΓ 先做 Γ，所以 AT 最后做 ΓT；Abc=PA+E 在 A 之后才用 P，所以转置后 P 要在 AT 之前用。前向是循环之后再覆盖：P(Av)+Ev；伴随必须循环之前清零、循环之后加回。照搬前向的做法，在循环之后令 ūi←wi，算出来的是 P(ATw)+Ew，它不是转置：下面的点积测试会失败，第一组随机向量的相对误差是 0.18。
Figure 9.1 sets the three functions side by side. The order is the whole point. A=Bu+BgΓ applies Γ first, so AT applies ΓT last; Abc=PA+E applies P after A, so after transposing P comes before AT. The forward code overwrites after the loop, P(Av)+Ev; the adjoint must zero before the loop and add back after it. Copying the forward recipe, setting ūi←wi after the loop, computes P(ATw)+Ew, which is not the transpose: the dot-product test below fails, with a relative error of 0.18 on the first random pair.
| residual | residual_d | residual_b | |||||
|---|---|---|---|---|---|---|---|
| 执行顺序 | order | 第一遍 → 第二遍 → 替换 | pass 1 → pass 2 → replace | 第一遍 → 第二遍 → 覆盖 | pass 1 → pass 2 → overwrite | 清零 → 第二遍T → 第一遍T → 加回 | zero → pass 2T → pass 1T → add back |
| 第一遍 | pass 1 | Gi = Σs wis(us − ui) | Ġi = Σs wis(u̇s − u̇i) | ū += ΓTḠ （读 Ḡ，写 ū）(reads Ḡ, writes ū) | |||
| 第二遍 | pass 2 | r = Σ ± h | ṙ = Σ ± ḣ | h̄ = r̄i − r̄j → ū, Ḡ | |||
| 壁面行 | wall rows | ri ← ui − uw | ṙi ← u̇i | 之前 r̄i ← 0；之后 ūi += wibefore: r̄i ← 0; after: ūi += wi |
图 9.2 逐句执行 residual_b。三行数组依次是输出 y（即 ū）、梯度的伴随 G_b（即 Ḡ）和输入 w̃（清零之后的 r̄）。第二遍倒着走时，每个面只把状态槽位直接写进 y，梯度槽位先存进 G_b；等所有面走完，第一遍倒着走，才把 G_b 经由权重 wis 送回 y。选壁面节点的 ê1 作种子，可以看到它在循环前被清零、整个循环什么也不写，最后只在 ū1 上加回 1：AbcT 的第 1 列只有对角线上的 1。
Figure 9.2 executes residual_b statement by statement. The three arrays are the output y (ū), the gradients' adjoint G_b (Ḡ) and the input w̃ (r̄ after the clearing). While pass 2 runs backwards, each face writes only its state slots into y directly and parks its gradient slots in G_b; only when every face is done does pass 1 run backwards and send G_b into y through the weights wis. Seed with the wall node's ê1 to see it zeroed before the loop, the whole loop write nothing, and a single 1 added back into ū1 at the end: column 1 of AbcT holds only the 1 on its diagonal.
上面把两遍循环写成矩阵再转置。同一件事也可以在代码的层面做：residual_d 的每一句都是几个切线之间的一个小线性映射，转置它、再把语句的顺序倒过来，就得到 residual_b。只需要三条规则（左边是前向的一句，右边是它在反向扫描里变成的那一句）：
Above, the two passes were written as matrices and transposed. The same can be done at the level of the code: every statement of residual_d is a small linear map between a few tangents; transpose each one and reverse the order of the statements, and residual_b comes out. Three rules are all it takes (left, a forward statement; right, what it becomes in the reverse sweep):
点是切线，横线是伴随（第 8 节的记号；代码里写作 t_d、t_b，而 u̇、ṙ、r̄、ū 是数组 v、r_d、w、y）。读变成累加（原来读 u̇，现在往 ū 里加），线性组合把伴随按同样的系数 η1、η2 分回各输入，累加写变成读（原来往 ṙ 里加，现在读 r̄）。每一条都是对偶恒等式在一句代码上的样子：输出的伴随乘输出的切线，等于各输入的伴随乘各自的切线之和。
A dot marks a tangent and a bar an adjoint (the notation of Section 8; the code writes t_d and t_b, and keeps u̇, ṙ, r̄, ū in the arrays v, r_d, w, y). A read becomes an accumulation (where u̇ was read, ū now accumulates), a linear combination hands its adjoint back to its inputs with the same coefficients η1, η2, and an accumulating write becomes a read (where ṙ was added to, r̄ is now read). Each is the duality identity written for one statement: the output's adjoint times its tangent equals the sum, over the inputs, of each input's adjoint times its tangent.
第一步，第二遍的最后两句（第 13、14 行）ṙi+=ḣ、ṙj−=ḣ，按规则 (c) 变成 residual_b 第 7 行（r̄ 已清零壁面分量，见第五步）：
Step 1, the last two statements of pass 2 (lines 13 and 14), ṙi+=ḣ and ṙj−=ḣ, become line 7 of residual_b by rule (c) (r̄ with its wall components cleared, step 5):
第二步，第 12 行 ḣ=(cL+νp)u̇L+(cR−νp)u̇R−νĠe·k，按规则 (b) 把 h̄ 分回三个输入；再把第 11 行 Ġe=12(Ġi+Ġj) 同样转置，就是 residual_b 第 8 至 10 行：
Step 2, line 12, ḣ=(cL+νp)u̇L+(cR−νp)u̇R−νĠe·k, hands h̄ back to its three inputs by rule (b); transposing line 11, Ġe=12(Ġi+Ġj), the same way gives lines 8 to 10 of residual_b:
状态直接累加进 ū（规则 (a)：原来读 u̇ 的地方现在往 ū 里累加）；梯度的份额先存进 Ḡ。
The states go straight into ū (rule (a): where u̇ was read, ū now accumulates); the gradients' share is parked in Ḡ.
第三步，第一遍（第 4 至 6 行）Δ̇=u̇j−u̇i、Ġi+=wijΔ̇、Ġj−=wjiΔ̇（Δ 是第 5 节第一遍跨过一个面的差值），倒着转置：先按 (b) 收集 Δ̄，再按 (a) 分给两个状态：
Step 3, pass 1 (lines 4 to 6), Δ̇=u̇j−u̇i, Ġi+=wijΔ̇ and Ġj−=wjiΔ̇ (Δ is the difference across a face of pass 1, Section 5), transposed backwards: first collect Δ̄ by rule (b), then hand it to the two states by rule (a):
这就是 residual_b 第 12 至 15 行的 γi−γj。
This is the γi−γj of lines 12 to 15 of residual_b.
第四步，顺序。前向先走完第一遍、再走第二遍；转置后语句的顺序倒过来，第二遍的转置要先对所有面走完，Ḡ 才完整，第一遍的转置才能开始。矩阵的写法里是同一件事（图 9.3）：
Step 4, the order. The forward code finishes pass 1 before pass 2; transposed, the order of the statements reverses: pass 2 transposed must run over every face first, so that Ḡ is complete before pass 1 transposed can start. In matrix form it is the same statement (Figure 9.3):
Bu 是第二遍对状态的部分，Bg 是对梯度的部分（上文的记号），Γ 是第一遍。Bu 正是一阶部分 A1（两者逐元素相同），梯度路径全在乘积 BgΓ 里。
Bu is the part of pass 2 acting on the states and Bg the part acting on the gradients (the symbols of the text above), Γ is pass 1. Bu is exactly the first-order part A1 (the two agree entry by entry), and the whole gradient path is in the product BgΓ.
第五步，壁面行。residual_d 最后的覆盖 ṙi←u̇i（第 16 行）其实是两句：先清零 ṙi←0（乘 P），再累加 ṙi+=u̇i（加 Eu̇）。它在前向里是最后一句，转置后就成了第一句：清零变成 r̄i←0，即“循环之前把 w 的壁面分量清零”（第 2 行），累加按规则 (a)、(c) 变成 ūi+=wi，即“把 wi 加回 ūi”（第 17 行；加回与循环一样只往 ū 里累加，先后无关，代码放在最后）。这正是 AbcTw=AT(Pw)+Ew（图 9.4）。
Step 5, the wall rows. The final overwrite of residual_d, ṙi←u̇i (line 16), is really two statements: clear ṙi←0 (multiply by P), then accumulate ṙi+=u̇i (add Eu̇). Last in the forward code, it comes first once transposed: the clearing becomes r̄i←0, “zero the wall components of w before the loops” (line 2), and the accumulation, by rules (a) and (c), becomes ūi+=wi, “add wi back into ūi” (line 17; like the loops it only adds into ū, so its place does not matter, and the code puts it last). That is exactly AbcTw=AT(Pw)+Ew (Figure 9.4).
核对：用构建时生成的一对伪随机向量做点积测试，⟨w,Abcv⟩ 与 ⟨AbcTw,v⟩ 的相对差是 3.6×10−16；residual_b 与显式的 AbcTw 最多差 1.3×10−15；按上式由三个因子重组的 A 与 jacobian() 最多差 8.9×10−16，Bu 与 A1 逐位相同（本页代码在构建时计算）。
Check: a dot-product test with a pair of pseudo-random vectors generated at build time gives a relative difference of 3.6×10−16 between ⟨w,Abcv⟩ and ⟨AbcTw,v⟩; residual_b differs from the explicit AbcTw by at most 1.3×10−15; the A reassembled from the three factors above differs from jacobian() by at most 8.9×10−16, and Bu agrees with A1 to the last bit (computed by this page's code at build time).
上面先按矩阵、再按三条规则把 residual_d 转置成了 residual_b。这里逐句做完：第 8 节每一句的切线都配上它的伴随，边界半面与设计输入也在内，并按反向扫描真正执行的顺序排好（图 9.5）。
Above, residual_d was transposed into residual_b first as matrices and then by three rules. Here it is done statement by statement: every statement's tangent of Section 8 gets its adjoint, boundary half-faces and design inputs included, in the order in which the reverse sweep actually runs (Figure 9.5).
第一步，一句的伴随。设一句算 y=φ(x1,…,xm)，它的切线是 ẏ=Σk∂φ∂xkẋk。它的伴随把 ȳ 按同样的偏导数分回每个输入：
Step 1, the adjoint of one statement. Let a statement compute y=φ(x1,…,xm); its tangent is ẏ=Σk∂φ∂xkẋk. Its adjoint hands ȳ back to every input with the same partial derivatives:
右边一式就是对偶恒等式在一句上的样子：输入各得的份额乘各自的切线，加起来等于输出的伴随乘输出的切线。为什么各句要倒序执行，为什么每个横线量都先清零、再用 += 累加，本节末尾的推导（“从切向量到伴随：转置把语句顺序反过来”）对一个面的四条语句讲清楚了。本页里，一个状态被它所有的面读，一个梯度被两侧的面读，读到它的每一句都要把自己的那一份加进它的伴随。
The right-hand formula is the duality identity for one statement: the shares the inputs receive, times their tangents, add up to the output's adjoint times its tangent. Why the statements run in reverse, and why every bar is cleared first and then accumulated with +=, is explained for the four statements of a face in the derivation at the end of this section (“From tangent to adjoint: the transpose runs the statements backwards”). On this page a state is read by all its faces and a gradient by the faces on both sides, and every statement that reads one adds its own share to its adjoint.
第二步，第二遍倒过来，逐面（residual_b 第 7 至 10 行）。最后执行的分发最先转置：ri+=h 与 rj−=h 变成读 r̄（壁面分量换成 −J̄，第四步）：
Step 2, pass 2 backwards, face by face (lines 7 to 10 of residual_b). The scatters ran last, so they are transposed first: ri+=h and rj−=h become reads of r̄ (the wall components replaced by −J̄, step 4):
再转置通量这一句：两部分各按自己的偏导数把 h̄ 分回输入，几何量也在内：
Then the flux statement: each part hands h̄ back to its inputs with its own partial derivatives, the geometric ones included:
最后转置收集：平均梯度把 Ḡe 两端各分一半，Ḡi+=12Ḡe、Ḡj+=12Ḡe；右状态的伴随回到它的来处：内部对偶面进 ūj，入流半面进 ū∞，出流半面进 ūi（uR=uL，所以出流半面上 ūi 共得 (cL+cR)h̄）。residual_b 只要 ū，把 ū∞、n̄、ξ̄ 这几份丢掉；第 11、12 节求设计导数时正是把它们收集起来。
Last, the gathers transposed: the averaged gradient gives each end half, Ḡi+=12Ḡe and Ḡj+=12Ḡe; the right state's adjoint goes back where the state came from: into ūj on an interior dual face, into ū∞ on an inflow half-face, into ūi on an outflow half-face (there uR=uL, so ūi receives (cL+cR)h̄ in all). residual_b wants ū only, and drops the shares of ū∞, n̄ and ξ̄; the design derivatives of Sections 11 and 12 collect exactly those.
第三步，第一遍倒过来，逐边（第 12 至 15 行）。两句累加 Gi+=wijΔ、Gj+=−wjiΔ 先收集 Δ̄，差值 Δ=uj−ui 再把它分给两个状态；权重的伴随留给第 12 节：
Step 3, pass 1 backwards, edge by edge (lines 12 to 15). The two accumulations Gi+=wijΔ and Gj+=−wjiΔ first collect Δ̄, and the difference Δ=uj−ui then hands it to the two states; the weights' adjoints are for Section 12:
第四步，源项与壁面行。Ri←Ri−si 转置成 s̄i−=r̄i。壁面行 ri←ui−uw,i 最后执行，所以最先转置：它的 r̄i 不经过通量（循环之前清零，第 2 行），直接给 ūi+=r̄i（第 17 行）与 ūw,i−=r̄i。种子取 r̄=ψ 时，后者就是第 11 节的 dJ/duw,i=−ψi。目标函数 J=−Σi∈WRi 读替换之前的壁面行，所以带着 J̄ 扫描时，循环看到的壁面分量是 −J̄：种子取 r̄=ψ、J̄=1 时，循环看到的正是第 12 节的 φ；源项的 s̄i−=r̄i 读的也是它，所以壁面节点上 s̄i=1。
Step 4, the source and the wall rows. Ri←Ri−si transposes into s̄i−=r̄i. The wall rows ri←ui−uw,i ran last and are transposed first: their r̄i does not pass through the fluxes (cleared before the loops, line 2) and goes straight into ūi+=r̄i (line 17) and ūw,i−=r̄i. Seeded with r̄=ψ, the latter is Section 11's dJ/duw,i=−ψi. The objective J=−Σi∈WRi reads the wall rows before the replacement, so a sweep that carries J̄ gives the loops −J̄ in the wall components: with r̄=ψ and J̄=1 the loops see exactly the φ of Section 12, and the source's s̄i−=r̄i reads the same value, so s̄i=1 on a wall node.
图 9.5 把三列并排：左列是 residual 的语句，中列是切线（从上往下执行），右列是伴随（从下往上执行：右列自下而上读，就是 residual_b 自上而下的伪代码）。图 9.6 取面 e20（第 3 节画过的那个面，这张网格上最歪的内部面之一），把这些语句在本页 u* 上的数值全部写出来。第 8 节末尾（“本页一个面的算例：前向扫描”，图 8.4）与本节末尾（“同一个面的反向扫描”，图 9.10）两处推导取的是另一个面 e17（节点 9 与 10 之间），种子相同，都是 v=du/du∞ 与 ψ。图 9.7 把 residual_b 整个跑一遍，种子取第 11 节的伴随解 ψ：ψ 满足 AbcTψ=−g，所以跑完以后 ū=−g，全部 24 个节点上都不为零。
Figure 9.5 sets the three columns side by side: on the left the statements of residual, in the middle their tangents (run top to bottom), on the right their adjoints (run bottom to top: read from the bottom up, the right column is the pseudocode of residual_b from the top down). Figure 9.6 takes the face e20 (the face drawn in Section 3, one of the most skewed interior faces of this mesh) and writes out every one of these statements' numbers at this page's u*. The derivations at the end of Section 8 (“A worked example on this page: the forward sweep of one face”, Figure 8.4) and at the end of this section (“The reverse sweep of the same face”, Figure 9.10) work another face, e17 (between nodes 9 and 10), with the same seeds, v=du/du∞ and ψ. Figure 9.7 runs residual_b once through, seeded with the adjoint solution ψ of Section 11: ψ satisfies AbcTψ=−g, so at the end ū=−g, non-zero on all 24 nodes.
核对：第一至四步逐句写成程序以后，用状态、来流值、源项、壁值、坐标上的随机种子与随机的 r̄、J̄，对偶恒等式 r̄Tṙ+J̄J̇=ūTu̇+D̄TḊ 两边的相对差最多 9.6×10−16；只对状态时与 residual_b 的相对差最多 2.4×10−16；种子取 r̄=ψ、J̄=1，一次扫描得到的 ū=AbcTψ+g 最大 3.6×10−16，D̄ 与本页代码分别算出的 dJ/du∞、dJ/ds、dJ/duw、dJ/dX 最多差 4.4×10−16；面 e20 上的对偶恒等式两边差 2.0×10−18；动画最后一帧与 −g 最多差 2.9×10−16（本页代码在构建时计算）。
Check: with steps 1 to 4 written as a program statement by statement, random seeds on the state, the free-stream value, the sources, the wall values and the coordinates, and random r̄ and J̄, the two sides of the duality identity r̄Tṙ+J̄J̇=ūTu̇+D̄TḊ differ by at most 9.6×10−16 relative; on the state alone the program differs from residual_b by at most 2.4×10−16 relative; seeded with r̄=ψ and J̄=1, one sweep gives ū=AbcTψ+g no larger than 3.6×10−16, and a D̄ within 4.4×10−16 of the dJ/du∞, dJ/ds, dJ/duw and dJ/dX this page's code computes separately; the two sides of the duality identity on the face e20 differ by 2.0×10−18; and the animation's last frame differs from −g by at most 2.9×10−16 (computed by this page's code at build time).
| 行 | line | residual 的语句 | the statement of residual | 切线 residual_d：从上往下执行 | tangent, residual_d: runs top to bottom | 伴随 residual_b：从下往上执行 | adjoint, residual_b: runs bottom to top |
|---|---|---|---|---|---|---|---|
| 第一遍：梯度，逐边 (i,j) | pass 1: the gradients, edge (i,j) by edge | ||||||
| 4 | Δ=uj−ui | Δ̇=u̇j−u̇i | ūj+=Δ̄ ūi−=Δ̄ （第 13、15 行）(lines 13, 15) | ||||
| 5 | Gi+=wijΔ | Ġi+=wijΔ̇+ẇijΔ | Δ̄+=wij·Ḡi w̄ij+=ΔḠi （第 12 行）(line 12) | ||||
| 6 | Gj+=−wjiΔ | Ġj+=−wjiΔ̇−ẇjiΔ | Δ̄−=wji·Ḡj w̄ji−=ΔḠj （第 14 行）(line 14) | ||||
| 第二遍：通量，逐个对偶面（壁面半面跳过） | pass 2: the fluxes, dual face by dual face (wall half-faces skipped) | ||||||
| 10 | uL=ui uR=uj 边界半面：uR 取 u∞ 或 uLboundary half-face: uR is u∞ or uL | u̇L=u̇i u̇R=u̇j | ūi+=ūL ūj+=ūR （第 8、9 行）(lines 8, 9) | ||||
| 11 | Ge=12(Gi+Gj) | Ġe=12(Ġi+Ġj) | Ḡi+=12Ḡe Ḡj+=12Ḡe （第 10 行）(line 10) | ||||
| 12 | ha=14vn(uL2+uR2)−ε|n|(uR−uL) | ḣa=cLu̇L+cRu̇R+∂ha∂n·ṅ | ūL+=cLh̄ ūR+=cRh̄ n̄+=∂ha∂nh̄ （第 8、9 行）(lines 8, 9) | ||||
| 12 | hv=−ν[p(uR−uL)+Ge·k] 远场半面：hv=0far field: hv=0 | ḣv=−ν[p(u̇R−u̇L)+Ġe·k]+∂hv∂n·ṅ+∂hv∂ξ·ξ̇ | ūL+=νph̄ ūR−=νph̄ Ḡe=−νkh̄ n̄,ξ̄+=… （第 8、9、10 行）(lines 8, 9, 10) | ||||
| 12 | h=ha+hv | ḣ=ḣa+ḣv | h̄a=h̄ h̄v=h̄ | ||||
| 13 | ri+=h | ṙi+=ḣ | h̄+=r̄i （第 7 行）(line 7) | ||||
| 14 | rj−=h | ṙj−=ḣ | h̄−=r̄j （第 7 行）(line 7) | ||||
| 循环之后：源项与壁面行 | after the loops: the source and the wall rows | ||||||
| ri−=si | ṙi−=ṡi | s̄i−=r̄i | |||||
| 16 | ri←ui−uw,i i∈W；伴随的循环在这里看到的是 −J̄（residual_b 里是 0）i∈W; the adjoint's loops see −J̄ here (0 in residual_b) | ṙi←u̇i−u̇w,i | ūi+=r̄i ūw,i−=r̄i （第 2、17 行）(lines 2, 17) | ||||
| 量 | quantity | 式子 | formula | 数值 | value |
|---|---|---|---|---|---|
| primal：本页的 u* | primal: this page's u* | ||||
| uL,uR | 0.1516, 0.1802 | ||||
| n | (0.2882, 0.3195) | ||||
| |n|,vn | vn=β·n | 0.4303, 0.4479 | |||
| ξ | xR−xL | (0.04630, 1.366) | |||
| p | ξ·n/|ξ|2 | 0.2409 | |||
| k | n−pξ | (0.2770, −0.009393) | |||
| GL,GR | (0.07305, −0.02100), (−0.05435, 0.001009) | ||||
| Ge | 12(GL+GR) | (0.009346, −0.009994) | |||
| cL,cR | 12vnuL+ε|n|,12vnuR−ε|n| | 0.3782, −0.3038 | |||
| ha,hv,h | −0.003653, −0.004793, −0.008446 | ||||
| 切线：种子是第 10 节的 v=du/du∞ | tangent: the seed is Section 10's v=du/du∞ | ||||
| u̇L,u̇R | 0.1853, 0.2343 | ||||
| ĠL,ĠR | (0.08339, −0.01526), (−0.07925, 0.01000) | ||||
| Ġe | 12(ĠL+ĠR) | (0.002073, −0.002630) | |||
| ḣa | cLu̇L+cRu̇R | −0.001139 | |||
| ḣv | −ν[p(u̇R−u̇L)+Ġe·k] | −0.006208 | |||
| ḣ | ḣa+ḣv | −0.007347 | |||
| 伴随：种子是第 11 节的 ψ | adjoint: the seed is Section 11's ψ | ||||
| r̄L,r̄R | ψL,ψR | −0.7804, −0.6758 | |||
| h̄ | r̄L−r̄R | −0.1046 | |||
| ūL | +=(cL+νp)h̄ | −0.05215 | |||
| ūR | +=(cR−νp)h̄ | 0.04437 | |||
| ḠL,ḠR | +=−12νkh̄ | (0.007243, −2.46×10⁻⁴) | |||
| h̄ḣ | 7.68×10⁻⁴ | ||||
| Σx̄ẋ | ūLu̇L+ūRu̇R+ḠL·ĠL+ḠR·ĠR | 7.68×10⁻⁴ | |||
| 几何：第 12 节的反向扫描，这个面的种子 h̄e | geometry: Section 12's reverse sweep, this face's seed h̄e | ||||
| h̄e | φi−φj | −0.1046 | |||
| ∂h∂n | ∂ha∂n+∂hv∂n | (−0.006682, −0.02041) | |||
| ∂hv∂ξ | (−0.001856, 0.002590) | ||||
| n̄ | h̄e∂h/∂n | (6.99×10⁻⁴, 0.002134) | |||
| ξ̄ | h̄e∂hv/∂ξ | (1.94×10⁻⁴, −2.71×10⁻⁴) | |||
对任意 v、w，恒等式 ⟨w,Abcv⟩=⟨AbcTw,v⟩ 必须在舍入误差内成立。用第 8 节的记号，取 u̇=v、r̄=w，于是 ṙ=Abcu̇、ū=AbcTr̄，这就是反向模式的对偶恒等式 r̄Tṙ=ūTu̇。图 9.8 用两组随机向量检验 residual_d 与 residual_b。第一组向量的相对误差是 1.9e−16；把壁面的转置做法换成“照搬前向”，同一组向量的误差变成 0.18。
For any v and w the identity ⟨w,Abcv⟩=⟨AbcTw,v⟩ must hold to round-off. In the notation of Section 8, with u̇=v and r̄=w, so that ṙ=Abcu̇ and ū=AbcTr̄, this is the duality identity of reverse mode, r̄Tṙ=ūTu̇. Figure 9.8 tests residual_d and residual_b on two random vectors. On the first pair the relative error is 1.9e−16; switch the wall's transpose to "copy the forward" and, on the same pair, it becomes 0.18.
第一步，测试的是什么。设 residual_b 算的是 Bw，B 是某个矩阵，写对了就是 AbcT。点积测试比较的两边之差是一个双线性型：
Step 1, what the test measures. Say residual_b computes Bw for some matrix B, which is AbcT if the code is right. The difference of the two sides the test compares is a bilinear form:
第二步，写对与写错。B=AbcT 时差只剩舍入误差。Abc−BT≠0 时，对固定的 w（只要 (Abc−BT)Tw≠0），wT(Abc−BT)v=0 只在一个超平面上成立，随机取的 v 落在上面的概率是零；差的大小与 Abc−BT 相对 Abc 的大小同一个量级。所以一对随机向量就足以判断：不用形成矩阵，也没有步长。
Step 2, right and wrong. With B=AbcT only round-off is left. With Abc−BT≠0, for a fixed w (as long as (Abc−BT)Tw≠0) the form wT(Abc−BT)v vanishes only on a hyperplane, which a random v hits with probability zero; its size is that of Abc−BT relative to Abc. So one random pair decides it: no matrix is formed, and there is no step size.
第三步，几种常见的写错，每种只改 residual_b 的一处（图 9.9）：
Step 3, some common slips, each changing one place of residual_b (Figure 9.9):
核对：每一列 30 对随机向量，residual_b 的各列用的是第 12 节表中那一组。本页的 residual_b 相对差最多 3.1×10−14（就是那张表里的数），几何的一对最多 3.1×10−14；几种写错里最小的也有 8.4×10−5，比舍入误差高出 9 个数量级以上（本页代码在构建时计算）。
Check: 30 random pairs per column, the residual_b columns on the pairs of the table of Section 12. This page's residual_b stays within 3.1×10−14 relative (the number in that table), and the geometric pair within 3.1×10−14; the smallest error among the slips is 8.4×10−5, more than 9 orders of magnitude above round-off (computed by this page's code at build time).
把所有残差和状态分量展平后，下面只是交换两个有限求和的次序。残差种子 是任意给定的，尚未求解伴随方程；只有选择满足第 11 节线性方程的 ，才得到该目标对应的离散伴随解。取 、，前向得到 ；以 为种子反向扫描，得到 。下式就是 AD 的对偶（点积）恒等式：
After flattening all residual and state components, the following simply exchanges two finite sums. The residual seed is arbitrary and has not solved an adjoint equation; only satisfying Section 11 gives the discrete adjoint solution for the chosen objective. With and the forward sweep gives ; the reverse sweep seeded with gives . The identity below is the duality (dot-product) identity of AD:
格点页的残差代码在面循环与源项之后还有一条语句（5）：把 8 个壁面节点的守恒行覆盖为 。 是保留行的对角掩码， 在被替换行的对角位置放 1；覆盖之前的残差与 Jacobian 就是正文的 与 ， 是它对设计变量的导数，面循环与源项写出的就是它们。
The residual code of the node page has one more statement (5) after the face loop and the sources: it overwrites the balance rows of the 8 wall nodes by . masks the retained rows and puts 1 on the diagonal of the replaced rows; the residual and Jacobian before the overwrite are the and of the main text, is its derivative with respect to the design, and the face loop and the sources write exactly these.
这条语句最后执行，所以它的反向语句（5′）最先执行：种子 分成两部分，保留行的部分 交给转置的面循环与源项，壁面行的部分 交给约束：对每个壁面节点 令 与 。完整转置必须包含约束的 Jacobian，而不是简单地把壁面伴随设为零；被替换的壁面行也不能把种子交给面循环与源项，所以壁面节点上的源项强度 只通过目标函数进入梯度。
This statement runs last, so its reverse statement (5′) runs first: the seed splits into the part of the retained rows, , which goes to the transposed face loop and the sources, and the part of the wall rows, , which goes to the constraints: for every wall node , and . The transpose includes the constraint Jacobian; it does not simply zero the wall adjoints, and the replaced wall rows must not hand their seeds to the face loop or to the sources, so the source strength of a wall node enters the gradient only through the objective.
内部面通量向左加、向右减，所以先把这两条路径合并成通量的伴随 ；之后按通量、重构、原始状态的相反依赖顺序传播。反向语句中的 += 很重要：同一状态会被多个面和多个求积点读取，这些贡献必须累加。这里没有分部积分，也没有另行猜测连续伴随边界条件。
An interior flux is added on the left and subtracted on the right, so combine those two paths into its adjoint before reversing fluxes, reconstruction and stored states. The += is essential: multiple faces and quadrature points read the same state, and their contributions must accumulate. There is no integration by parts and no separately guessed continuous adjoint boundary condition.
第 8 节局部链里的每个切向量都有对应的伴随。反向语句按与前向相反的顺序执行，每个输入的横线量都用 += 累加：
Every tangent of the local chain of Section 8 has its adjoint. The reverse statements run in the opposite order of the forward ones, and the bar of every input accumulates with +=:
一个面的切向量是四个初等 Jacobian 的乘积。乘积的转置等于各个因子的转置按相反的顺序相乘：
The tangent of one face is the product of four elementary Jacobians. The transpose of a product is the product of the transposes in the reverse order:
仍从右往左读，这就是反向扫描的语句顺序：程序最后执行的分发最先处理，最先执行的取值最后处理。先把所有伴随置零，再以 为种子逐条执行；覆盖（5）在整个面循环之后执行，所以它的反向语句（5′）排在最前：
Read again from right to left, this is the order of the reverse sweep: the scatter, which the program ran last, comes first, and the gather, which it ran first, comes last. Zero every adjoint, seed with , and run; the overwrite (5) ran after the whole face loop, so its reverse statement (5′) comes first:
面循环结束后，源项的一对也按同样的规则反过来：，，；被替换的行已由上面的（5′）处理：它们的种子不进入面循环，写进这些行的源项因此得不到权重。本页代码里每一条语句的伴随语句（顺序与前向相反）也都在本节里逐句写出。
After the face loop the source is reversed by the same rule: , , ; the replaced rows were handled by (5′) above: their seeds do not enter the face loop, so a source written into such a row receives no weight. The adjoint statement of every statement of this page's code, in the reverse order, is written out in this section as well.
边界面执行同样的语句。远场入流面上 的右槽位不是存储状态，而是 ，所以 的这一行为零：（2′）中 不把它的伴随 传给 ， 则把它加进 的 分量。
A boundary face runs the same statements. On a far-field inflow face the right slot of is not a stored state but , so its row of is zero: in (2′), passes none of its adjoint to , and adds it to the component of .
分发与取值互为转置。分发 把一个面通量写进两行；它的转置 从这两行读回一个数，是一次取值。取值 从全局数组读出几个状态；它的转置 把它们的伴随加回原来的位置，是一次带累加的分发。
Scatter and gather are transposes of each other. The scatter writes one face flux into two rows; its transpose reads those two rows back into one number: a gather. The gather reads a few states out of the global array; its transpose adds their adjoints back into their places: a scatter with accumulation.
为什么要先置零、再用 += 累加：程序里的一次扇出（一个值被多处读取）在反向扫描里变成扇入。一个存储状态属于它每个相邻面的 ，设计变量也进入许多面与源项；链式法则对变量的每一次使用各给出一项， 是所有这些项的和：
Why zero first and then +=: a fan-out of the program (one value read in several places) becomes a fan-in of the reverse sweep. A stored state belongs to the of every neighbouring face, and the design enters many faces and the sources; the chain rule gives one term per use of a variable, and is the sum of all these terms:
让这个和从零开始，每个 += 加上一次使用的贡献。反过来，程序里的扇入（ 收集它各个面的通量）在反向扫描里变成扇出： 只被读取，不需要累加。
starts this sum at zero, and each += adds the term of one use. Conversely, a fan-in of the program ( collecting the fluxes of all its faces) becomes a fan-out of the reverse sweep: is only read and needs no accumulation.
把同一个面的前向语句与反向语句配对。每一步只用一次转置的定义 ，把一个初等 Jacobian 从内积的右边移到左边，所以内积在每一级都相同：
Pair the forward and reverse statements of the same face. Each step uses the definition of the transpose once, , to move one elementary Jacobian from the right of the inner product to the left, so the inner product is the same at every stage:
对所有面求和，再加上源项的一对 （），右边的两项正好累加成反向扫描的输出 与 ：
Summed over all faces, with the pair of the source added ( with ), the two terms on the right accumulate exactly into the outputs and of the reverse sweep:
被替换的行再加上它们约束的一对 （ 是切向量一节里覆盖语句的切向量）：它正是（5′）里约束的状态路径与设计路径对 的贡献；所以同一个恒等式对本页的方程组也成立（ 与 取整个反向扫描结束时的值），下面的测试用的就是它：
The replaced rows add the pair of their constraints, ( is the tangent of the overwrite in the tangent section): it is exactly what the state path and the design path of the constraints in (5′) contribute to . The same identity therefore holds for this page's system ( and taken at the end of the whole reverse sweep), and the test below uses that system:
这就是 AD 的点积测试：任取种子 与 ，前向扫描算出 ，反向扫描算出 ，两边之差应在舍入误差量级；只要有一条反向语句的转置写错，差就是 的。本页在构建时用它自己的 matrix-free 前向与反向乘积（ 个分量的确定性种子，设计方向取下面算例的 ）做了这个测试：，，相对差 。
This is the dot-product test of AD: take any seeds and , run the forward sweep for and the reverse sweep for ; the two sides must agree to round-off, and a single wrong transpose in the reverse statements gives an difference. This page ran the test on its own matrix-free forward and reverse products at build time (deterministic seeds with components, the design direction the of the example below): and , a relative difference of .
反向种子取本页伴随求解的解 ；面循环的种子是掩码后的 ，这个面的两行都是保留行，所以它在这两行等于 。于是 与 是这个面对 与 的贡献（也是对 与 的贡献）。下表按反向语句的顺序列出伴随值；最后一列是每一级的内积：先是 ，然后是 、，最后是 ，切向量取第 8 节“本页一个面的算例：前向扫描”里同一个面的前向算例的值。
Seed the reverse sweep with this page's adjoint solution, ; the seed of the face loop is the masked one, , and both rows of this face are retained, so it equals in them. Then and are this face's contributions to and (and so to and ). The table lists the adjoints in the order of the reverse statements; the last column is the inner product at each stage: first , then , , and last , with the tangents of the forward example of the same face, the one under “A worked example on this page: the forward sweep of one face” in Section 8.
| 语句Statement | 量Quantity | 伴随Adjoint | 内积Inner product |
|---|---|---|---|
| （0）置零(0) zero | 0 | ||
| （4′）种子(4′) seed | −0.7804 | −1.573816158×10⁻⁴ | |
| −0.7942 | |||
| （4′）取值(4′) gather | 0.01384 | −1.573816158×10⁻⁴ | |
| （3′）通量的转置(3′) flux, transposed | 0.008203 | −1.573816158×10⁻⁴ | |
| −0.008768 | |||
| 8.863×10⁻⁴ 3.556×10⁻⁴ | |||
| 8.863×10⁻⁴ 3.556×10⁻⁴ | |||
| （2′）（1′）加回(2′) (1′) add back | −4.82×10⁻⁴ | −1.573816158×10⁻⁴ | |
| −4.583×10⁻⁴ | |||
| 0.008473 | |||
| −0.00828 | |||
| −5.141×10⁻⁴ | |||
| −4.587×10⁻⁴ | |||
| 5.891×10⁻⁴ | |||
| 7.066×10⁻⁴ | |||
| −1.397×10⁻⁴ | |||
| 0 |
4 级内积相同，最大相对差 ：对偶恒等式在这些数上成立；图 9.10 按反向语句逐条播放这些数。这个面不直接依赖 ，所以 ； 的伴随来自它真正进入的那些行。
The 4 inner products agree, with a largest relative difference of : the duality identity holds in these numbers. Figure 9.10 plays them statement by statement in reverse. This face does not depend on directly, so ; the adjoint of comes from the rows it does enter.
箭头全部反向，每个初等 Jacobian 转置，顺序相反；存储状态与设计变量处用 += 累加。下方：转置的乘积与形状。
Every arrow reversed, every elementary Jacobian transposed, in the reverse order; the stored states and the design accumulate with +=. Below: the transposed product and its shapes.
（0）置零并给种子：，； −0.7804， −0.7942（这个面的两行都是保留行，掩码后的种子 在这两行等于 ）。
(0) Zero and seed: , ; −0.7804, −0.7942 (both rows of this face are retained, so the masked seed equals in them).
（4′）分发的转置是一次取值： 0.01384。内积 −1.573816158×10⁻⁴。
(4′) The transposed scatter is a gather: 0.01384. Inner product −1.573816158×10⁻⁴.
（3′）通量的转置：， 0.008203， −0.008768， (8.863×10⁻⁴, 3.556×10⁻⁴)， (8.863×10⁻⁴, 3.556×10⁻⁴)。内积 −1.573816158×10⁻⁴。
(3′) The transposed flux: , 0.008203, −0.008768, (8.863×10⁻⁴, 3.556×10⁻⁴), (8.863×10⁻⁴, 3.556×10⁻⁴). Inner product −1.573816158×10⁻⁴.
（2′）重构的转置：，，这个面给出 0（）。
(2′) The transposed reconstruction: and ; this face gives 0 ().
（1′）取值的转置是带累加的分发：， −4.82×10⁻⁴， −4.583×10⁻⁴， 0.008473， −0.00828， −5.141×10⁻⁴， −4.587×10⁻⁴， 5.891×10⁻⁴， 7.066×10⁻⁴， −1.397×10⁻⁴。内积 −1.573816158×10⁻⁴。
(1′) The transposed gather is a scatter with accumulation: , −4.82×10⁻⁴, −4.583×10⁻⁴, 0.008473, −0.00828, −5.141×10⁻⁴, −4.587×10⁻⁴, 5.891×10⁻⁴, 7.066×10⁻⁴, −1.397×10⁻⁴. Inner product −1.573816158×10⁻⁴.
对偶：4 级内积相同（最大相对差 5×10⁻¹⁵）， 在这个面上逐级成立。
Duality: the 4 inner products agree (largest relative difference 5×10⁻¹⁵); holds stage by stage on this face.
第 8 节只算了乘积 Abcv。前向模式真正要解的是切线方程：选设计变量的一个分量 Dk（种子 Ḋ=ek），对收敛解满足的 Rbc=0 求导，u 随 Dk 的变化率 v=du/dDk，也就是这个种子下状态的切线 u̇，满足
Section 8 only computed the product Abcv. What forward mode actually has to solve is the tangent equation: pick one design component Dk (the seed Ḋ=ek) and differentiate the converged residual equation Rbc=0 with respect to it; the rate v=du/dDk at which u follows Dk, which is the state's tangent u̇ for this seed, satisfies
这里取远场来流值 Dk=u∞。只有远场入流半面见到它（作为对流通量的右状态；远场没有扩散通量），所以 f=∂Rbc/∂u∞ 就是这些半面上通量对右槽位的导数，写进各自的节点，共 5 个非零分量。壁面行 ui−uw 与 u∞ 无关，那里 fi=0。
Here take the far-field free-stream value, Dk=u∞. Only the far-field inflow half-faces see it, as the right state of their convective flux (the far field carries no diffusive flux), so f=∂Rbc/∂u∞ is the derivative of those fluxes with respect to their right slot, written into their nodes: 5 non-zero entries. The wall rows ui−uw do not involve u∞, so fi=0 there.
切线方程是线性的，用第 6 节那个 Jacobi 迭代、同一个 diag(Abc) 求解：
The tangent equation is linear, and it is solved by the Jacobi iteration of Section 6, with the same diag(Abc):
每一轮只调用一次与矩阵有关的函数：第 8 节的 residual_d，算出 Abcv。diag(Abc) 是第 6 节的 diagonal 在 u* 处的值。除此之外，矩阵既不形成，也不分解：这就是 matrix-free 的前向求解。壁面行上 (Abc)ii=1、(Abcv)i=vi、fi=0，更新就是 vi←vi−vi=0：壁值不随 u∞ 变，v 在壁面上每一轮都是 0。
Each sweep calls just one matrix-related function: residual_d of Section 8, which computes Abcv. diag(Abc) is what diagonal of Section 6 returns at u*. Beyond that the matrix is neither formed nor factorised: this is the forward solve, matrix-free. On a wall row (Abc)ii=1, (Abcv)i=vi and fi=0, so the update reads vi←vi−vi=0: the wall values do not move with u∞, and v is 0 on the wall in every sweep.
在解附近，primal 迭代的误差每一轮乘以 I−diag(Abc)−1Abc，切线迭代的误差乘的正是同一个矩阵。所以 primal 一过它非线性的起始阶段，两条收敛曲线就以同一个速率下降，这个速率就是谱半径 0.787：实测 k=60 到 100 之间，两者每轮都缩小 0.7874 倍。起点不同：primal 从 u≡1 出发，切线从 v=0 出发，初始残差就是 f。图 10.1 逐轮播放切线迭代，看切线解 v 如何从远场入流半面铺进区域。切线迭代在 k=121 时降到判据 10−13 以下，从 k=155 起落在舍入误差平台上。
Near the solution the primal iteration multiplies its error every sweep by I−diag(Abc)−1Abc, and the tangent iteration multiplies its error by exactly the same matrix. So once the primal is past its nonlinear start, the two convergence curves fall at the same rate, the spectral radius 0.787: measured between k=60 and 100, both shrink by 0.7874 per sweep. They start from different places: the primal from u≡1, the tangent from v=0, its first residual being f. Figure 10.1 plays the tangent iteration sweep by sweep, showing how the tangent solution v spreads into the domain from the far-field inflow half-faces. The tangent iteration drops below the 10−13 criterion at k=121 and is on its round-off floor from k=155 on.
入流半面上对流通量的右状态就是 u∞，所以 f 只来自这一项，写进它左侧的节点：
On an inflow half-face the right state of the convective flux is u∞ itself, so f comes from that term alone, written into the node on its left:
入流半面上 vn<0，两项都是负的：f 在 5 个节点上为负，其余为零。解 Av=−f（这里是 Abc）得到的 v 在非壁面节点上处处为正（图 10.2）：来流值抬高，整个区域的 u 都跟着抬高，只有壁面节点不动。
On an inflow half-face vn<0, so both terms are negative: f is negative in 5 nodes and zero elsewhere. Solving Av=−f (here with Abc) gives a v that is positive at every node off the wall (Figure 10.2): raise the free-stream value and u rises over the whole domain, except on the wall nodes, which stay put.
核对：把 u∞ 挪动 ±10−6 重算残差（状态不动），中心差分与上式的 f 最多差 8.1×10−11（本页代码在构建时计算）。
Check: moving u∞ by ±10−6 and re-evaluating the residual (state held fixed), central differences differ from the f above by at most 8.1×10−11 (computed by this page's code at build time).
J 不直接依赖 u∞，所以
J does not depend on u∞ directly, so
本算例是 3.69011506499406；在 u∞±10−6 各重解一次 primal 的中心差分与它的相对差是 2.2×10−10。一次切线求解只回答一个设计分量：要得到对第 12 节 48 个坐标的导数，前向那一路要解 48 次。下一节从另一头走，一次伴随求解就得到全部。
Here it comes to 3.69011506499406; a central difference that re-solves the primal at u∞±10−6 differs from it by 2.2×10−10 (relative). One tangent solve answers one design component: for the derivatives with respect to the 48 coordinates of Section 12, the forward route has to solve 48 times. The next section goes from the other end, where one adjoint solve gives them all.
第一步，上面解的 u∞ 是第 1 节设计变量 D 的 81 个分量之一。每个分量 Dk 的切线方程与上面的相同，只是右端项 fk 不同，目标函数还带着自己的直接项：
Step 1, the u∞ solved above is one of the 81 components of the design variables D of Section 1. Every component Dk has the tangent equation above, with its own right-hand side fk, and the objective adds its own direct term:
第二步，右端项从哪几句来，第 8 节第六步已经列出。按类读出来：
Step 2, which statements each right-hand side comes from is listed in step 6 of Section 8. Read by kind:
第三步，都用同一个 Jacobi 迭代、同一个 diag(Abc) 求解。迭代矩阵相同，收敛速率也相同，只是起点的残差 fk 不同：到判据要 84 到 121 轮（壁面节点上的源项 f=0，一轮不用）。图 10.3 每类取一个，画出 f 与解出的 v。
Step 3, all are solved by the same Jacobi iteration with the same diag(Abc). The iteration matrix is the same, and so is the rate; only the starting residual fk differs: 84 to 121 sweeps to the criterion (a source on a wall node has f=0 and needs none). Figure 10.3 takes one of each kind and draws f and the solution v.
核对：81 个分量各解一次，gTvk+∂J/∂Dk 与第 12 节一次伴随扫描给出的 D̄k 最多差 3.3×10−16（最大的分量是 3.690；本页代码在构建时计算）。
Check: with each of the 81 components solved once, gTvk+∂J/∂Dk and the D̄k of one adjoint sweep (Section 12) differ by at most 3.3×10−16 (the largest component is 3.690; computed by this page's code at build time).
设原始残差已经收敛，且当前光滑分支上的 非奇异。取设计种子（方向），令 ；前向求解得到的状态切向量是 。对每一个残差分量逐项求导，零残差的切向量仍为零：
Assume the primal residual is converged and is nonsingular on the current smooth branch. Choose a design seed (direction) and set ; the state tangent returned by the tangent solve is . Differentiate each residual component: the tangent of a zero residual is zero.
实现分三步：状态种子 、设计种子 ，前向求值残差得到 ；用精确 的乘积求解 ；再以 和 为种子对目标前向求导，得到 。一次 只是算子调用（一次前向扫描），不是求解 。稳态导数的对象是收敛残差方程，不是 Jacobi 的迭代次数。在不可微点（如 limiter 的切换或分支判别处）以上普通导数需要另外讨论，不能无条件套用。
Implementation has three steps: seed the state with and the design with and evaluate the residual tangent ; solve for with exact products; then differentiate the objective forward with the seeds and to obtain . One is an operator application (one forward sweep), not the solution . The steady derivative concerns the converged residual equation, not the Jacobi iteration count. At nondifferentiable points (a limiter switch or a branch test) these ordinary derivatives need separate treatment.
用面循环的语句来看：右端 是 时的一次前向扫描。它在保留的行上是面循环的 ，只有 与源项两条路径起作用；在被替换的行上是约束自己的设计路径，即壁面行里的 ，任何面循环里都没有它：对壁面值 ，整个右端都来自这条路径（覆盖之前的残差不依赖它们）。迭代中的每个 是 时的一次前向扫描：保留的行是面循环的，被替换的行是约束对状态的路径。一次前向求解的全部工作就是这段面循环与覆盖以不同的种子反复执行。对 个设计分量逐个求导，要做 次这样的求解。
In terms of the face loop's statements: the right-hand side is one forward sweep with . In the retained rows it is the face loop's , in which only the paths and the source are active; in the replaced rows it is the constraint's own design path, in the wall rows, which no face loop contains: for the wall values the whole right-hand side comes from that path (the residual before the overwrite does not depend on them). Every of the iteration is one forward sweep with : the face loop in the retained rows and the constraint's state path in the replaced rows. All the work of a tangent solve is this face loop and overwrite run again and again with different seeds, and the derivatives with respect to all design components take such solves.
伴随方程 AbcTψ=−g 的右端项 g=(∂J/∂u)T 由目标函数决定（第 7 节）。它同样用第 6 节的 Jacobi 迭代、同一个 diag(Abc) 求解（转置不改变对角元）：
The right-hand side of the adjoint equation AbcTψ=−g, g=(∂J/∂u)T, is set by the objective (Section 7). It too is solved by the Jacobi iteration of Section 6, with the same diag(Abc) (transposition leaves the diagonal alone):
与矩阵有关的唯一操作是一次第 9 节的 residual_b 调用：它的种子是 r̄=ψ，收敛时返回 ū=AbcTψ=−g。替换行也自动对上了：Abc 的壁面行是 êiT，(Abc)ii=1，转置之后是第 i 列上的 êi，正是 residual_b 里“清零、加回”所对应的结构。
The only matrix-related operation is one call to residual_b of Section 9: its seed is r̄=ψ, and at convergence it returns ū=AbcTψ=−g. The replaced rows match up by themselves: the wall rows of Abc are êiT with (Abc)ii=1, which transposed become êi in column i, exactly the structure the zeroing and adding back of residual_b implement.
伴随解同样给出第 10 节的那个导数，不用再解任何方程：
The adjoint solution gives the derivative of Section 10 as well, with no further solve:
本算例两者都是 3.69011506499406，相对差 1.2×10−16。一次伴随求解回答一个输出对所有设计分量的导数：同一个 ψ 也给出第 12 节对 48 个坐标的导数，那里前向那一路要解 48 次。
Here both come to 3.69011506499406, with a relative difference of 1.2×10−16. One adjoint solve answers the derivatives of one output with respect to every design component: the same ψ also gives the derivatives with respect to the 48 coordinates of Section 12, where the forward route solves 48 times.
伴随迭代的误差每一轮乘以 I−diag(Abc)−1AbcT，它是 I−Abcdiag(Abc)−1 的转置，而
The adjoint iteration multiplies its error every sweep by I−diag(Abc)−1AbcT, the transpose of I−Abcdiag(Abc)−1, and
与 primal、切线的 diag(Abc)−1Abc 相似。转置不改变特征值，所以三个迭代矩阵的谱完全相同：进入渐近阶段以后，primal、切线、伴随三条收敛曲线以同一个速率下降，这个速率就是谱半径 0.787。实测 k=60 到 100 之间，三条曲线每轮都缩小 0.7874 倍，四位有效数字相同。三条曲线的起点不同：primal 从 u≡1 出发，切线和伴随从 0 出发，初始残差就是 f 与 g；也只有 primal 有一段非线性的起始阶段。但一进入渐近阶段，它们就平行。图 11.1 把三条曲线画在一起，并逐轮播放伴随迭代，看伴随解 ψ 如何从壁面铺开：有扩散，上游、下游都有。
is similar to the primal's and the tangent's diag(Abc)−1Abc. Transposition leaves eigenvalues unchanged, so the three iteration matrices have exactly the same spectrum: once in the asymptotic phase the primal, tangent and adjoint convergence curves fall at the same rate, the spectral radius 0.787. Measured between k=60 and 100, all three shrink by 0.7874 per sweep, the same to four digits. The curves start from different places: the primal from u≡1, the tangent and the adjoint from 0, so that their first residuals are f and g; and only the primal has an early nonlinear phase. But once asymptotic they run parallel. Figure 11.1 draws the three curves together and plays the adjoint iteration sweep by sweep, showing how the adjoint solution ψ spreads out from the wall: with diffusion, downstream as well as upstream.
第一步，为什么 ψTf=gTv。前向方程给出 f=−Abcv，伴随方程给出 AbcTψ=−g；把前者代入、再用转置的定义挪一下括号：
Step 1, why ψTf=gTv. The tangent equation gives f=−Abcv and the adjoint equation AbcTψ=−g; substitute the first and move the brackets by the definition of the transpose:
这是代数恒等式，与网格粗细无关；两个数的差只来自两个迭代各自的舍入误差。
It is an algebraic identity, whatever the mesh; the two numbers differ only by the round-off of the two iterations.
第二步，前向迭代的误差。精确解 v* 满足 Abcv*+f=0，所以 v*=v*−diag(Abc)−1(Abcv*+f)；与迭代式相减：
Step 2, the error of the tangent iteration. The exact solution v* satisfies Abcv*+f=0, so v*=v*−diag(Abc)−1(Abcv*+f); subtract this from the iteration:
方程是线性的，这一式是精确的，不需要线性化：前向迭代的误差每一轮乘的正是 primal 渐近段的那个矩阵。
The equation is linear, so this is exact, with no linearisation: every sweep multiplies the tangent's error by exactly the matrix of the primal's asymptotic phase.
第三步，伴随迭代的误差。同样地，
Step 3, the error of the adjoint iteration. In the same way,
第四步，这个矩阵与前者的转置相似。diag(Abc) 是对角阵，(diag(Abc)−1)T=diag(Abc)−1，于是
Step 4, this matrix is similar to the transpose of the former. diag(Abc) is diagonal, (diag(Abc)−1)T=diag(Abc)−1, and so
相似变换不改变特征值，转置也不改变特征值（矩阵与它的转置有相同的特征多项式）。所以三个迭代矩阵的谱完全相同，谱半径相同：渐近段三条收敛曲线平行。它们的起点不同、初期不同，只有渐近段的斜率是同一个。
A similarity transformation leaves the eigenvalues alone, and so does transposition (a matrix and its transpose have the same characteristic polynomial). The three iteration matrices therefore have exactly the same spectrum and the same spectral radius: in the asymptotic phase the three convergence curves run parallel. They start from different places and differ early on; only the asymptotic slope is shared.
核对：本页 ψTf 与 gTv 的相对差是 1.2×10−16；第四步的恒等式两边逐元素最多差 5.6×10−17，两个矩阵的 24 个特征值逐一配对最多差 1.2×10−15；第 60 到 100 轮实测，primal、前向、伴随每轮分别乘 0.78741、0.78741、0.78741，谱半径 ρ=0.78741（本页代码在构建时计算）。图 11.2 把三条曲线画到舍入误差平台。
Check: on this page ψTf and gTv differ by 1.2×10−16 relative; the two sides of the identity of step 4 differ entry by entry by at most 5.6×10−17, and the 24 eigenvalues of the two matrices, paired one by one, by at most 1.2×10−15; measured over sweeps 60 to 100, the primal, the tangent and the adjoint shrink per sweep by 0.78741, 0.78741 and 0.78741, against the spectral radius ρ=0.78741 (computed by this page's code at build time). Figure 11.2 draws the three curves down to their round-off floors.
收敛到机器精度是什么样子。三个迭代都没有停在判据处，而是跑满 300 轮。primal、切线、伴随分别在 k=121、121、122 时降到判据 10−13 以下，但这时仍以每轮 0.787 倍的比例下降：判据只是约定的停止点，并不说明迭代已经走到头。它们按同一比例一直降到 10−16 附近，然后分别从 k=162、155、156 起落在舍入误差平台上：此后再也不超过这次运行最后三分之一里的最大值，三个平台都在 7.0×10−17 与 1.9×10−16 之间。曲线变平，才说明迭代已经收敛到机器精度。停在判据处的 ψ 与平台上的 ψ 相差 2.5×10−13。
What convergence to machine precision looks like. None of the three iterations stops at the criterion; each runs the full 300 sweeps. The primal, the tangent and the adjoint drop below the 10−13 criterion at k=121, 121 and 122, but there they are still falling by a factor of 0.787 per sweep: a criterion is an agreed place to stop, not a sign that the iteration has run its course. They keep falling at that rate down to about 10−16 and then settle onto their round-off floors from k=162, 155 and 156 on: from there they never again rise above the largest value of the run's last third, and the three floors all lie between 7.0×10−17 and 1.9×10−16. Only this flat tail shows that an iteration has converged to machine precision. The ψ at the criterion is 2.5×10−13 away from the ψ on the floor.
令 φ 等于 ψ，只是在壁面节点上取 −1（图 11.3）。第 12 节还会用到它。它满足替换之前的伴随方程：对每一个不在壁面上的节点 j，(ATφ)j=0（实测最大 4.1e−16）。这是 g=−Σi∈WAi,:T 的直接结果：伴随方程在非壁面节点上是 Σk∉WAkjψk=Σi∈WAij，把右边移过去，就是 φ 在壁面上取 −1 的齐次方程。
Let φ equal ψ except on the wall nodes, where it is −1 (Figure 11.3); Section 12 needs it too. It satisfies the adjoint equation before the replacement: (ATφ)j=0 for every node j off the wall (measured: at most 4.1e−16). This follows directly from g=−Σi∈WAi,:T: off the wall the adjoint equation reads Σk∉WAkjψk=Σi∈WAij, and moving the right-hand side over gives the homogeneous equation for φ with −1 on the wall.
它的物理意义很直接。在节点 i 的控制体里加一个源项 si，则 dJ/dsi=−φi：−φi 就是这个单位源最终流进物面的比例，其余的从远场流走。壁面节点上是 1（全部）；中圈节点在 0.71 到 0.79 之间；外圈节点在 0.49 到 0.68 之间。对流页面上，伴随场主要集中在物体上游：不在壁面上的下游节点只有 0.01 到 0.07 左右，那是耗散项带过去的；这里有了物理扩散，φ 在所有非壁面节点上都在 −0.79 与 −0.49 之间，下游也不例外。
Its physical meaning is direct. Add a source si in the control volume of node i; then dJ/dsi=−φi, so −φi is the fraction of that unit source that ends up in the body, the rest leaving through the far field. It is 1 on the wall nodes (all of it), between 0.71 and 0.79 on the middle ring, and between 0.49 and 0.68 on the outer ring. On the advection page the adjoint field sits mostly upstream of the body: the downstream nodes off the wall carry only about 0.01 to 0.07, brought there by the dissipation term; here, with physical diffusion, φ lies between −0.79 and −0.49 on every node off the wall, downstream ones included.
第一步，壁面以外的节点。强度为 si 的源项进入节点 i 的平衡：Ri←Ri−si（第 15 节的源项，一点求积）。它是设计变量 D 的一个分量；替换后的 Rbc 同样看到它，∂Rbc/∂si=−êi，而 J=−Σk∈WRk 不含它。于是由伴随公式：
Step 1, a node off the wall. A source of strength si enters the balance of node i as Ri←Ri−si (the source term of Section 15, one-point quadrature). It is one component of the design variables D; the replaced Rbc sees it the same way, ∂Rbc/∂si=−êi, and J=−Σk∈WRk does not contain it. The adjoint formula gives
第二步，壁面节点。源项只进入替换之前的壁面行 Ri；求解用的行已换成 ui−uw，所以状态根本不动。但目标函数 J=−Σk∈WRk 直接含有这一行：
Step 2, a wall node. The source enters only the wall row Ri before the replacement; the solve uses ui−uw in its place, so the state does not move at all. But the objective J=−Σk∈WRk contains that row directly:
这就是 φ 在壁面上取 −1 的原因：壁面节点上的源项全部被物面吸收。
This is why φ takes −1 on the wall: a source on a wall node is absorbed by the body entirely.
第三步，为什么是“份额”。把所有节点的平衡加起来：内部对偶面的通量两两抵消，剩下物面、远场与源项。记 J∞ 为从远场流出的净通量，则在收敛解上
Step 3, why it is a share. Add up the balances of all nodes: the fluxes through interior dual faces cancel in pairs, and what is left is the body, the far field and the sources. With J∞ the net flux out through the far field, at the converged solution
一个单位源分成两份：−φi 进物面，1+φi 从远场出去。
A unit source splits in two: −φi into the body and 1+φi out through the far field.
核对：在 n1、 n9、 n17 各放一个 ±10−6 的源项并重解原始方程，dJ/ds 依次是 1.0000、0.7804、0.6376，与 −φ 最多差 4.6×10−10；从远场出去的份额依次是 0.0000、0.2196、0.3624，两份之和与 1 最多差 5.3×10−10（本页代码在构建时计算）。
Check: a source of ±10−6 in n1, n9 and n17 in turn, re-solving the primal each time, gives dJ/ds= 1.0000, 0.7804 and 0.6376, within 4.6×10−10 of −φ; the shares that leave through the far field are 0.0000, 0.2196 and 0.3624, and the two shares add up to 1 within 5.3×10−10 (computed by this page's code at build time).
壁面节点上的 ψi 本身也有意义。把壁值看作参数，替换行 ui−uw,i 对它的偏导是 −1，所以
On the wall nodes ψi has a meaning of its own. Take the wall values as parameters: the replaced row ui−uw,i has derivative −1 with respect to its own, so
本算例中 ψi 在 0.197 与 0.303 之间：把某个壁面节点的壁值提高一点，流进物面的通量就相应减少，这符合直觉。它与中心差分（每次重解原始方程）的最大差是 5.0e−10。同一次伴随求解同时给出了对形状（第 12 节）和对壁值的灵敏度。
In this test case ψi lies between 0.197 and 0.303: raising the wall value at one wall node lowers the flux into the body accordingly, as expected. It agrees with central differences (re-solving the primal each time) to 5.0e−10. One adjoint solve gives the sensitivities to the shape (Section 12) and to the wall values together.
−φi 是节点 i 处一个单位源流进物面的份额（上文）。壁面节点上它是 1：源就在物面上，全部被吸收。这正是附录 A.1 第 3 小节连续伴随的壁面条件 λ=1（附录 A.1 第 1 小节定义 λh=−ψh，附录 A.3 用它与离散伴随比较，壁面节点上取 −φi=1）。离开物面，份额下降多快，取决于伴随怎样输运。
−φi is the share of a unit source at node i that ends up in the body (above). On a wall node it is 1: a source on the body is absorbed entirely. That is the wall condition λ=1 of the continuous adjoint of Appendix A.1, part 3 (Appendix A.1, part 1 defines λh=−ψh, and Appendix A.3 compares the discrete adjoint with it through that relation, taking −φi=1 on the wall nodes). How fast the share falls away from the body depends on how the adjoint is transported.
伴随逆着来流输运。伴随方程的矩阵是 AbcT。在 Abc 里，一个面的对流通量主要读迎风一侧的状态，一个状态的变化顺着 β 传到下游的行；转置把这种依赖倒过来，伴随在一点上的值听的是它下游的伴随，于是伴随逆着 β 传。连续问题里也是这样：线性化的对流 uβ·∇(·) 转置以后是 −∇·(uβ·)，方向相反。而 u 在物面附近很小（物面上是 0），所以贴壁的一层里伴随主要靠扩散，离开物面、u 变大以后，逆流的输运才显出来。primal 正好相反：u 顺着 β 从远场带进来，迎流一侧在物面前形成一层薄的边界层，背流一侧留下长长的尾迹。
The adjoint is transported against the flow. The matrix of the adjoint equation is AbcT. In Abc a face's convective flux reads mostly the upwind state, and a change of a state travels along β into the rows downstream; transposing reverses that dependence, so the adjoint at a point listens to the adjoint downstream of it, and it travels against β. The continuous problem does the same: the linearised convection uβ·∇(·) transposes into −∇·(uβ·), the other way round. And u is small near the body (0 on it), so within the layer next to the wall the adjoint moves mostly by diffusion; only further out, where u has grown, does the upstream transport show. The primal is the mirror image: u is carried in from the far field along β, with a thin boundary layer in front of the body on the side facing the flow and a long wake behind it.
图 11.4 在第 15 节的加密网格族上解同一个伴随方程（第 0 至 3 层，每层都用本页的 Jacobi 迭代解到舍入误差平台），沿两条从物体中心出发的射线取值。离物面 0.5 处，上游对下游的 −φ（第 0 至 3 层）依次是 0.742 对 0.780、0.771 对 0.766、0.831 对 0.749、0.840 对 0.740：从第 1 层起上游更大，即伴随逆着来流伸得更远；第 0 层反过来。u 则是每一层物面的影响都顺着来流伸得更远：第 3 层上游已恢复到 0.622，下游只有 0.120。
Figure 11.4 solves the same adjoint equation on the mesh family of Section 15 (levels 0 to 3, each by this page's Jacobi iteration down to its round-off floor) and samples it along two rays from the body's centre. At 0.5 from the body, −φ upstream against downstream (levels 0 to 3) is 0.742 against 0.780, 0.771 against 0.766, 0.831 against 0.749 and 0.840 against 0.740: from level 1 on the upstream value is larger, so the adjoint reaches further against the flow; on level 0 it is the other way round. For u the body's influence reaches further along the flow on every level: on level 3 it has recovered to 0.622 upstream and is only 0.120 downstream.
核对：每一层存下的 ψ 代回伴随方程，AbcTψ+g 的均方根最多 1.0×10−16；ψTf 与同一层切线解的 gTv（附录 A.3 表中的 dJ/du∞）相对差最多 4.0×10−16（本页代码在构建时计算）。
Check: put back into the adjoint equation, each level's stored ψ leaves AbcTψ+g with an rms of at most 1.0×10−16; and ψTf differs from the same level's tangent gTv (the dJ/du∞ of the table in Appendix A.3) by at most 4.0×10−16 relative (computed by this page's code at build time).
在满足 的解上 。先把 固定，沿切向量 、 对 前向求导，再按 和 收集系数。即使 随设计改变，额外的 项也在收敛解上为零。
At a feasible state , . Hold fixed, differentiate forward along the tangents and , then collect the coefficients of and . Even if varies with the design, the extra term vanishes at the converged state.
选择 使未知的状态敏感度前面的系数为零。于是一个目标只需求解一次伴随线性系统，再进行设计反传；总设计梯度是列向量 。注意符号来自这里的加号约定，并不是因为“伴随一定带负号”。
Choose to cancel the coefficient of the unknown state sensitivity. One objective then requires one adjoint solve followed by design reverse accumulation; the total design gradient is a column. Its sign follows from the plus convention in this Lagrangian, not from a universal negative sign for adjoints.
用反向模式的横线记号：以 为种子运行残差的反向扫描，得到 ，这就是伴随方程。再把目标自己的反向语句（种子 ）累加进同一组横线量：它们给 加上 ，给 加上 ，于是 。伴随方程说的正是状态的横线量总和为零：状态不再把敏感度往下传，同一次扫描的设计伴随 就是梯度 。也可以直接检查与 tangent 的等价性：利用 ，，所以两条路线给出同一个 。更换为 会得到 ，但这只是两个离散乘子的符号约定，不能据此认定连续场 与离散向量逐点相等。
In the bar notation of reverse mode: the residual's reverse sweep seeded with returns , which is the adjoint equation. Accumulate the objective's own reverse statements (seeded with ) into the same bars: they add to and to , so . The adjoint equation says exactly that the state's total bar vanishes: the state passes no sensitivity on, and the design adjoint of the same sweep, , is the gradient . Check equivalence with the tangent: gives , so both routes give the same . Changing to yields , but that relates two discrete sign conventions; it does not prove pointwise equality with a continuous field .
精确算子决定求解哪个方程；Jacobi 迭代决定如何求解。较便宜的近似 Jacobian（例如一阶或冻结部分系数的 Jacobian）可以用于迭代修正或预处理，但不能替代精确转置。应监控真残差 。若其不为零，当前伴随与精确 tangent 的方向导数差为 ，因此只看迭代次数或预处理残差不够。
The exact operator determines the equation; the Jacobi iteration determines how it is solved. A cheaper approximate Jacobian (first order, or with some coefficients frozen) may supply corrections or preconditioning, but cannot replace the exact transpose. Monitor the true residual . If it is nonzero, the current adjoint directional derivative minus the exact tangent derivative is ; an iteration count or a preconditioned residual alone is insufficient.
同样的结论也可以不经过 Lagrangian，直接由前向求解的两个关系得到：
The same result follows without a Lagrangian, directly from the two relations of the tangent solve:
最后一项是三个因子 、、 的乘积，括号可以从两边加。前向路线先算 ：它依赖方向 ，每个设计分量都要解一次以 为系数的方程，全部分量要解 次。伴随路线先算行向量 ：它与 无关，只要解一次。把这个行向量记作 ：
The last term is the product of three factors, , and , which can be bracketed from either side. The tangent route computes first: it depends on the direction , so every design component needs one solve with , solves for all of them. The adjoint route computes the row first: it does not depend on , so it is solved once. Call this row :
这对每个方向 都成立，所以梯度就是上面方框里的 ：一次伴随求解之后，再以 为种子做一次反向扫描，就得到全部分量。两条路线的计算量（本页的 与 由本页代码数出，与页首路线图相同）：
This holds for every direction , so the gradient is the boxed above: after one adjoint solve, one reverse sweep seeded with gives every component. The cost of the two routes ( and of this page, counted by its code, as in the page roadmap at the top):
| 前向（tangent）路线Tangent route | 伴随（adjoint）路线Adjoint route | |
|---|---|---|
| 每个目标的线性求解Linear solves per objective | 次（其中 8 次右端为零），系数 (8 of them with a zero right-hand side), with | 1 次，系数 1, with |
| 右端Right-hand sides | ||
| 求解之后After the solves | 次前向扫描，各得一个 forward sweeps, one each | 1 次反向扫描，得到整个 1 reverse sweep, the whole |
| 一次求解给出One solve gives | 一个分量 one component | 全部 个分量all components |
| 个目标时With objectives | 仍是 次（ 与目标无关）still ( does not depend on the objective) | 次 |
本页 、：前向路线要解 81 个 24 阶的线性方程组（其中 8 个的右端 为零：壁面节点上的源项强度 ，它们所在的行被约束替换，所以实际要做的求解是 73 个），伴随路线只解 1 个。设计变量多而目标少时（形状优化通常如此）伴随路线便宜得多；设计变量少而目标多时，前向路线更合适。
This page has and : the tangent route solves 81 linear systems of order 24 (8 of them have a zero right-hand side : the source strengths of the wall nodes, whose rows are replaced by the constraint, so 73 solves do work), the adjoint route one. With many design variables and few objectives (the usual case in shape optimisation) the adjoint route is far cheaper; with few design variables and many objectives the tangent route is the better choice.
全导数 dJ/dX=∂J/∂X+ψT∂Rbc/∂X 里有两处要看清。第一，替换行 ui−uw 与坐标无关，它们在 ∂Rbc/∂X 里是零行，所以只有 ψ 的非壁面分量起作用。第二，J=−Σi∈WRi(u,X) 通过替换之前的壁面行依赖坐标，贡献正好是壁面上的 −1。两者合起来，权重就是第 11 节的 φ：
Two things in the total derivative dJ/dX=∂J/∂X+ψT∂Rbc/∂X need a close look. First, the replaced rows ui−uw do not depend on the coordinates: they are zero rows of ∂Rbc/∂X, so only the off-wall components of ψ act. Second, J=−Σi∈WRi(u,X) depends on the coordinates through the wall rows before the replacement, and contributes exactly −1 on the wall. Together, the weights are the φ of Section 11:
用反向模式的记号，这就是坐标的伴随 X̄=(∂Rbc/∂X)Tψ+(∂J/∂X)T，所以几何梯度就是替换之前的残差的一次反向扫描：残差的种子是 r̄=φ，每个面通量的伴随是 h̄e=φi−φj，把它经由几何量一路传回坐标。状态 u 固定在 u*。
In reverse-mode notation this is the adjoint of the coordinates, X̄=(∂Rbc/∂X)Tψ+(∂J/∂X)T, so the geometric gradient is one reverse sweep of the residual before the replacement: seeded with r̄=φ, each face flux has the adjoint h̄e=φi−φj, carried back through the geometric quantities to the coordinates, with the state held at u*.
一个面的通量通过四条路径依赖坐标（图 12.1）：对偶面法向 n（由边中点与单元中心构成，与对流页面相同）；边向量 ξ=xj−xi；由 n 与 ξ 得到的 p、k；以及经由模板向量 ξs 进入的最小二乘权重 wis。用 n=pξ+k 把扩散通量改写成 hv=−ν[Ge·n+p(Δu−Ge·ξ)]，它对 n 与 ξ 的偏导是：
A face's flux depends on the coordinates along four paths (Figure 12.1): the dual-face normal n (built from edge midpoints and cell centres, as on the advection page); the edge vector ξ=xj−xi; p and k, which come from n and ξ; and the least-squares weights wis, through the stencil vectors ξs. With n=pξ+k the diffusive flux reads hv=−ν[Ge·n+p(Δu−Ge·ξ)], and its partial derivatives with respect to n and ξ are
对流部分只经由 n：∂ha/∂n=14(uL2+uR2)β−ε(uR−uL)n/|n|，与对流页面相同。平均梯度 Ge 那一路给出梯度的伴随 Ḡi（与第 9 节的 Ḡ 同一个公式，Ḡi 是各面 −12νkh̄e 之和），再经由最小二乘传给模板向量。Gi=Mi−1Σsξs(us−ui)，Mi=ΣsξsξsT；对它求切线再转置，令 b̄=Mi−1Ḡi，得到
The convective part depends on n only: ∂ha/∂n=14(uL2+uR2)β−ε(uR−uL)n/|n|, as on the advection page. The averaged-gradient path yields the gradients' adjoints Ḡi (the formula of Section 9's Ḡ: Ḡi sums −12νkh̄e over the faces), which least squares passes on to the stencil vectors. With Gi=Mi−1Σsξs(us−ui) and Mi=ΣsξsξsT, take the tangent, transpose it and set b̄=Mi−1Ḡi:
推导：状态固定，对模板向量求切线，Ġi=Mi−1(Σsξ̇s(us−ui)−ṀiGi)，Ṁi=Σs(ξ̇sξsT+ξsξ̇sT)；代入对偶恒等式 Ḡi·Ġi=Σsξ̄s·ξ̇s 的左边，用 Mi 的对称性把 Mi−1 移到 Ḡi 上，再按 ξ̇s 整理即得。最后，ξ 与 ξs 都是坐标之差：x̄j+=ξ̄、x̄i−=ξ̄；对偶面法向的反向与对流页面第 12 节相同。
Derivation: with the state held fixed, the tangent with respect to the stencil vectors is Ġi=Mi−1(Σsξ̇s(us−ui)−ṀiGi) with Ṁi=Σs(ξ̇sξsT+ξsξ̇sT); substitute it into the left side of the duality identity Ḡi·Ġi=Σsξ̄s·ξ̇s, move Mi−1 onto Ḡi by the symmetry of Mi, and collect the terms by ξ̇s. Finally, ξ and ξs are both coordinate differences, x̄j+=ξ̄ and x̄i−=ξ̄; the reverse of the dual-face normals is as in Section 12 of the advection page.
一个坐标分量 Xm 并不牵动所有面。在 图 12.2 里点选一个节点，可以看到它移动时跟着变的几何量：法向 n 变化的对偶面（紫；对偶面由边中点与单元中心构成），边向量 ξ 从而 p、k 变化的面（虚线），最小二乘权重变化的节点（着色），以及通量因此变化的面（红）。反向扫描求 dJ/dXm 时，只有这些红色的面有贡献，每个面带着自己的面通量伴随 h̄e=φi−φj（标在面上）；紫色箭头就是这样得到的 dJ/dXm。
One coordinate Xm does not reach every face. Pick a node in Figure 12.2 to see what moves with it: the dual faces whose normal n moves (violet; a dual face is built from edge midpoints and cell centres), the faces whose edge vector ξ, and so p and k, move (dashed), the nodes whose least-squares weights move (tinted), and the faces whose flux moves as a result (red). In the reverse sweep for dJ/dXm only those red faces contribute, each with its own face adjoint h̄e=φi−φj (written on it); the violet arrow is the dJ/dXm they give.
第一步，依赖什么。Gi=Σswis(us−ui) 对状态是线性的；权重 wis=Mi−1ξs 只依赖几何（不加距离权，与 u 无关），经模板向量 ξs 依赖坐标。所以对状态求切线时权重不动，对坐标求切线时才要权重的切线。
Step 1, what depends on what. Gi=Σswis(us−ui) is linear in the state; the weights wis=Mi−1ξs depend on the geometry only (no distance weighting, no u), on the coordinates through the stencil vectors ξs. So a state tangent leaves the weights alone, and only a coordinate tangent needs theirs.
第二步，权重的切线。由 MiMi−1=I 得 ddtMi−1=−Mi−1ṀiMi−1，于是
Step 2, the weights' tangent. From MiMi−1=I, ddtMi−1=−Mi−1ṀiMi−1, and so
第三步，梯度的切线有两项：
Step 3, the gradient's tangent has two terms:
第一项是第 8 节的第一遍；第二项只在坐标带切线时出现。其中 ḃi=Σsξ̇s(us−ui) 是 bi=Σsξs(us−ui) 在状态不动时的切线，只含几何的切线；代入第二步，第二项就是上文状态固定时的 Ġi。
The first term is pass 1 of Section 8; the second appears only for a coordinate seed. In it ḃi=Σsξ̇s(us−ui) is the tangent of bi=Σsξs(us−ui) with the state held fixed, the geometric tangents only; with step 2, the second term is the Ġi at a fixed state above.
第四步，对线性场精确，网格移动时也不丢。对 ΣswisξsT=I 求切线，得 Σs(ẇisξsT+wisξ̇sT)=0。线性场 u=u0+∇u·x 在移动的点上取值时 us−ui=ξs·∇u、u̇s−u̇i=ξ̇s·∇u，代入第三步：
Step 4, exactness on linear fields survives a moving mesh. The tangent of ΣswisξsT=I is Σs(ẇisξsT+wisξ̇sT)=0. A linear field u=u0+∇u·x sampled at moving points has us−ui=ξs·∇u and u̇s−u̇i=ξ̇s·∇u, and step 3 gives
网格怎么动，线性场的梯度都还是 ∇u。
However the mesh moves, a linear field's gradient stays ∇u.
第五步，转置：梯度的伴随经同一组权重回到模板的各点。把第三步配成对偶恒等式 Ḡi·Ġi=Σsūsu̇s+ūiu̇i+Σsw̄is·ẇis，按各个切线读出系数：
Step 5, the transpose: the gradient's adjoint returns to the stencil's points through the same weights. Pair step 3 with the duality identity Ḡi·Ġi=Σsūsu̇s+ūiu̇i+Σsw̄is·ẇis and read off the coefficients of the tangents:
前两式是“转置的模板”：前向时 Gi 从模板各点收集，伴随时 Ḡi 按同样的权重分回这些点，每点一份，自己减去总和。这就是 residual_b 第 12 至 15 行的 γ。
The first two are the transposed stencil: forward, Gi gathers from the stencil's points; in the adjoint, Ḡi hands back to the same points with the same weights, one share each, and the node itself gives up their sum. These are the γ of lines 12 to 15 of residual_b.
第六步，再把第二步转置，从权重的伴随得到模板向量的伴随：
Step 6, transpose step 2, and the weights' adjoints give the stencil vectors':
代入第五步的 w̄is=Ḡi(us−ui)，并记 b̄=Mi−1Ḡi：Σsw̄isξsT=ḠibiT，所以 M̄i=−b̄GiT，三项合起来正是上文的 ξ̄s=b̄(us−ui)−b̄(Gi·ξs)−Gi(b̄·ξs)。
Put in step 5's w̄is=Ḡi(us−ui) and write b̄=Mi−1Ḡi: then Σsw̄isξsT=ḠibiT, so M̄i=−b̄GiT, and the three terms together are exactly the ξ̄s=b̄(us−ui)−b̄(Gi·ξs)−Gi(b̄·ξs) above.
第七步，一个节点的数值：节点 n9（中间一圈，5 个边邻居），Ḡi 取本节这一次扫描中收集到的 (0.0115,−0.0110)，b̄=Mi−1Ḡi=(0.0132,−0.0037)：
Step 7, the numbers of one node: node n9 (on the middle ring, 5 edge neighbours), with the Ḡi=(0.0115,−0.0110) collected by this section's sweep and b̄=Mi−1Ḡi=(0.0132,−0.0037):
| 模板点 s | stencil point s | ξs | wis | us−ui | wis·Ḡi | ξ̄s |
|---|---|---|---|---|---|---|
| n1n1 | (−0.4932, −0.01800) | (−0.5490, 0.01265) | −0.1516 | −0.006468 | (−0.001064, 2.96×10⁻⁴) | |
| n10n10 | (−0.3557, 0.8865) | (−0.4288, 0.2788) | −7.01×10⁻⁴ | −0.008018 | (0.001167, −3.32×10⁻⁴) | |
| n16n16 | (−0.4070, −0.8352) | (−0.4230, −0.2355) | 0.1256 | −0.002280 | (0.001993, −5.62×10⁻⁴) | |
| n17n17 | (0.6028, 0.02200) | (0.6709, −0.01547) | 0.06354 | 0.007905 | (−3.13×10⁻⁴, 9.15×10⁻⁵) | |
| n18n18 | (0.04630, 1.366) | (0.001624, 0.4077) | 0.02866 | −0.004477 | (0.001042, −2.95×10⁻⁴) |
核对：权重的切线与第 14 节的复数步长最多差 1.3×10−15；线性场放在随机移动的网格上，梯度的切线最大 4.8×10−16（第四步）；第五、六步包含在第 9 节的对偶恒等式核对里（本页代码在构建时计算）。
Check: the weights' tangent differs from the complex step of Section 14 by at most 1.3×10−15; for a linear field on a randomly moving mesh the gradient's tangent is at most 4.8×10−16 (step 4); steps 5 and 6 are part of the duality check of Section 9 (computed by this page's code at build time).
图 12.3 把这一次反向扫描逐步放出来。先逐面：每个面把自己的伴随 h̄e 经 ∂h/∂n、∂hv/∂ξ 变成 n 与 ξ 的伴随，推动构成面的各点和 ξ 的两端（对偶面由边中点与相邻单元的中心构成，这些点的推动再平均分给所在边或单元的顶点），梯度的份额存进 Ḡ。所有面走完、Ḡ 齐了以后，再逐节点：Ḡi 经最小二乘权重变成模板向量的伴随（上一段第六步），推动模板的各点。代码先把各点的伴随收集起来、最后再分给节点；动画每一步直接分，和是一样的。每一步的推动都成对出现、大小相等方向相反（n、ξ、ξs 都是两点之差），所以任何时刻全部箭头之和都停在舍入误差上：“整张网格平移不改变任何东西”逐面成立。最后一帧就是下面 图 12.4 的箭头。
Figure 12.3 plays this reverse sweep step by step. Face by face first: each face turns its adjoint h̄e, through ∂h/∂n and ∂hv/∂ξ, into adjoints of n and ξ that push the points the face is built from and the two ends of ξ (a dual face is built from an edge midpoint and the centres of the neighbouring cells, and each of these points shares its push equally among the nodes of its edge or cell), and parks its gradient share in Ḡ. Once every face is done and Ḡ is complete, node by node: Ḡi becomes the stencil vectors' adjoints through the least-squares weights (step 6 above) and pushes the stencil's points. The code collects the points' adjoints first and hands them to the nodes at the end; the animation hands them on at once, which gives the same sums. Every step pushes in pairs, equal and opposite (n, ξ and ξs are all differences of two points), so at any moment the sum of all arrows is at round-off: "translating the whole mesh changes nothing" holds face by face. The last frame is the arrows of Figure 12.4 below.
核对：动画最后一帧与本页 gradient() 的结果最多差 5.6×10−17；每一步推动的和最大 6.9×10−18，全部箭头之和在任何一帧最大 8.3×10−17（本页代码在构建时计算）。
Check: the animation's last frame differs from what this page's gradient() returns by at most 5.6×10−17; the pushes of one step sum to at most 6.9×10−18, and the sum of all arrows in any frame to at most 8.3×10−17 (computed by this page's code at build time).
图 12.4 画出这一次反向扫描的结果：每个网格节点上一个箭头 dJ/dXi。物面节点上的 8 个箭头全部朝外（与径向最多差 8.8°）：把物面往外推，物体吸收的通量变大。
Figure 12.4 draws the result of this reverse sweep: one arrow dJ/dXi at each mesh node. On the body nodes all 8 arrows point outwards (within 8.8° of the radial direction): push the body surface out, and the flux the body absorbs grows.
所有几何量都只依赖坐标之差，整张网格平移不改变任何东西，所以 ΣidJ/dXi 在两个方向上都必须为零。实测为 0 与 1.7e−16。转动则不然：对流速度 β 是固定的，转动网格会改变解。
Every geometric quantity depends on coordinate differences only, so translating the whole mesh changes nothing, and ΣidJ/dXi must vanish in both directions. Measured: 0 and 1.7e−16. Rotation is another matter: the convection velocity β is fixed, so rotating the mesh changes the solution.
| 精确 Jacobian 对中心差分（相对） | exact Jacobian against central differences (relative) | 3.0e−11 |
|---|---|---|
| residual_d、residual_b 对显式 Abc（第 9 节的那一对随机向量） | residual_d, residual_b against the explicit Abc (the random pair of Section 9) | 4.4e−16, 1.3e−15 |
| 点积测试（30 组随机向量中最差） | dot-product test (worst of 30 random pairs) | 3.1e−14 |
| (ATφ)j，j 不在壁面上 | (ATφ)j, j off the wall | 4.1e−16 |
| Jacobi 满足判据的轮次（primal / 切线 / 伴随） | Jacobi sweeps to the criterion (primal / tangent / adjoint) | 121 / 121 / 122 |
| dJ/du∞：切线 gTv 对伴随 ψTf（相对） | dJ/du∞: tangent gTv against adjoint ψTf (relative) | 1.2e−16 |
| dJ/dX 对全链路中心差分（h = 10⁻⁶，每次重解） | dJ/dX against full-chain central differences (h = 10⁻⁶, re-solving each time) | 绝对 1.1e−9，相对 3.4e−9absolute 1.1e−9, relative 3.4e−9 |
| Σi dJ/dXi | Σi dJ/dXi | (0, 1.7e−16) |
| dJ/duw,i = −ψi 对中心差分 | dJ/duw,i = −ψi against central differences | 5.0e−10 |
第一步，前向那一路：每个分量 Dk 一个切线方程（第 10 节末尾），共 81 个：
Step 1, the tangent route: one tangent equation per component Dk (the end of Section 10), 81 in all:
第二步，伴随那一路：一个伴随方程，再一次反向扫描。第 9 节逐句写出的反向扫描，种子取 r̄=ψ、J̄=1，一次就给出设计的伴随：
Step 2, the adjoint route: one adjoint equation and one reverse sweep. The reverse sweep written out statement by statement in Section 9, seeded with r̄=ψ and J̄=1, gives the design's adjoint in one go:
它的四部分都在前面：ū∞=ψTf（第 11 节），s̄i=−φi（第 11 节），壁值的 ūw（第 11 节），坐标的 X̄（本节）。
Its four parts all appear above: ū∞=ψTf (Section 11), s̄i=−φi (Section 11), the wall values' ūw (Section 11), and the coordinates' X̄ (this section).
第三步，两条路为什么给出同一个数。对第 k 个分量，把第一步代进第二步：
Step 3, why the two routes give the same number. For component k, put step 1 into step 2:
与第 11 节 u∞ 那一个分量的推导完全一样，只是对 81 个分量同时成立。
It is the derivation of Section 11 for the one component u∞, now for all 81 at once.
第四步，代价。与矩阵有关的计算都在 Jacobi 迭代里，每轮一次矩阵乘积。前向：81 个分量各有一个切线方程，各自解到判据，共 7754 轮（每轮一次 residual_d）；其中壁面节点上的 8 个源项 f=0，解就是 v=0，一轮不用，真正要解的是 73 个；伴随：一次求解 122 轮（每轮一次 residual_b），再加一次反向扫描。求解次数之比是 73 比 1（把那 8 个空的也算上是 81 比 1），轮数之比约 64 比 1。图 12.5 把两条路并排：伴随一步就立起全部 81 根柱，前向一个一个地把绿点落在同一根柱上，其中 8 个空的落下去不占轮数。
Step 4, the cost. Everything that involves the matrix sits in a Jacobi iteration, one matrix product per sweep. Tangent: 81 components, one tangent equation each, each solved to the criterion, 7754 sweeps in all, one residual_d each; the 8 sources on wall nodes have f=0, so their solution is v=0 and needs no sweep, and 73 of the solves are real ones; adjoint: one solve of 122 sweeps, one residual_b each, and one reverse sweep. The solves are 73 against 1 (81 against 1 if the 8 empty ones are counted), the sweeps about 64 against 1. Figure 12.5 sets the two routes side by side: the adjoint raises all 81 bars in one step, the tangent drops the green dots onto the same bars one at a time, the 8 empty ones without a sweep.
核对：两组 81 个数最多差 3.3×10−16（最大的分量是 3.690；本页代码在构建时计算）。
Check: the two sets of 81 numbers differ by at most 3.3×10−16 (the largest component is 3.690; computed by this page's code at build time).
上一节把 dJ/dX 推了出来：伴随只解一次线性方程，就拿到全部 48 个分量。这一节只做一件事：验证它是对的。判据是最朴素的那个，逐个坐标分量做中心差分：
The previous section derived dJ/dX: the adjoint solves one linear equation and gets all 48 components. This section does one thing only: check that it is right. The criterion is the plainest one there is, a central difference taken coordinate by coordinate:
关键在于 J(Xm±h) 是完整的目标函数：每次都要改动网格，重建全部几何量（对偶面法向、边向量、p、k、最小二乘权重），重新解一次非线性的原始方程，再算反力通量。所以这里验证的是全导数，状态对几何的响应也算在内。代价也正因如此：中心差分每个分量要两次非线性求解，一共 96 次，而伴随只用了一次线性求解。只往一侧扰动的前向差分把基准值 J(X) 共用，一共 49 次，代价是截断误差从 O(h2) 变成 O(h)。
What matters is that J(Xm±h) is the complete objective: every evaluation moves the mesh, rebuilds all the geometry (dual-face normals, edge vectors, p, k, least-squares weights), solves the nonlinear primal again, and only then computes the reaction flux. So what is verified is the total derivative, the response of the state to the geometry included. That is also what it costs: two nonlinear solves per component, 96 in all, against the single linear solve of the adjoint. A forward difference, perturbing to one side only, shares the baseline J(X) and needs 49, at the price of an O(h) truncation error instead of O(h2).
在 图 13.1 里可以自己改步长 h 和差分格式，再按“用差分验证”重算；底部逐个变量扫描步长。
In Figure 13.1 the step h and the scheme can be set by hand and recomputed with "Verify by differences"; the sweep at the foot runs through the step sizes one variable at a time.
最优步长在哪里。把 48 个分量逐一扫描：中心差分的最优步长有 40 个落在 h=10−5，其余在 10−4 或 10−6；前向差分有 44 个落在 10−7，其余在 10−6 或 10−8。这与 ϵmach1/3、ϵmach1/2 的估计（ϵmach 为机器精度）一致。各分量在最优步长下能达到的相对误差，中心差分的中位数是 6.0e−10，最好 1.2e−11；前向差分的中位数是 1.0e−7。默认的 h=10−6 对中心差分略偏小，对前向差分略偏大：在这个步长下，48 个分量中差分与伴随的最大绝对差，中心差分是 1.1e−9，前向差分是 1.6e−7。点积测试（第 9 节）则完全不含步长，是更干净的检查。差分的步长两难，在第 14 节的复数步长那里不存在。
Where the optimal step lies. Sweep the 48 components one by one: the central optimum is h=10−5 for 40 of them and 10−4 or 10−6 for the rest; the forward optimum is 10−7 for 44 and 10−6 or 10−8 for the rest. That matches the ϵmach1/3 and ϵmach1/2 estimates (ϵmach the machine precision). The relative error each component reaches at its optimal step has a median of 6.0e−10 for central differences, 1.2e−11 at best, and 1.0e−7 for forward differences. The default h=10−6 is a little small for the central scheme and a little large for the forward one: at that step the largest absolute difference from the adjoint over the 48 components is 1.1e−9 for central differences and 1.6e−7 for forward ones. The dot-product test (Section 9) has no step in it at all and is the cleaner check. The step dilemma of differences does not arise for the complex step of Section 14.
第一步，截断误差。记 J'、J''、J''' 为 J 对坐标 Xm 的各阶导数，由 Taylor 展开：
Step 1, the truncation error. With J', J'' and J''' the derivatives of J with respect to the coordinate Xm, Taylor expansion gives
中心差分里偶数阶项正负抵消，所以截断误差是 O(h2)；前向差分是 O(h)。
The even terms cancel in the central difference, so its truncation error is O(h2); the forward difference's is O(h).
第二步，舍入误差。每次算出的 J 都带着约 ϵJ≈ϵmach|J| 的误差（原始方程解到了舍入误差平台）。两个几乎相等的数相减，这个误差不会变小，再除以 2h 或 h：
Step 2, the round-off error. Every computed J carries an error of about ϵJ≈ϵmach|J| (the primal is solved to its round-off floor). Subtracting two nearly equal numbers does not shrink that error, and it is then divided by 2h or h:
第三步，求最小。令导数为零：
Step 3, minimise. Set the derivatives to zero:
代回去，最小误差分别正比于 ϵmach2/3 与 ϵmach1/2。取 ϵmach=2−52，ϵmach1/3=6.1×10−6，ϵmach1/2=1.5×10−8：步长比这更小，舍入误差占上风；更大，截断误差占上风。复数步长（第 14 节）没有相减，也就没有第二项。
Put back, the smallest errors are proportional to ϵmach2/3 and ϵmach1/2. With ϵmach=2−52, ϵmach1/3=6.1×10−6 and ϵmach1/2=1.5×10−8: at smaller steps round-off wins, at larger ones truncation. The complex step (Section 14) subtracts nothing, and so has no second term.
核对：图 13.2 取 dJ/dX 中最大的分量（节点 20 的 y 坐标），步长从 10−1 到 10−12 每次重解原始方程：中心差分在 h=10−4 处相对误差最小，为 3.3×10−11；前向差分在 h=10−7 处最小，为 3.0×10−8（本页代码在构建时计算）。
Check: Figure 13.2 takes the largest component of dJ/dX (the y coordinate of node 20) and re-solves the primal at every step from 10−1 to 10−12: central differences are best at h=10−4, with a relative error of 3.3×10−11; forward differences at h=10−7, with 3.0×10−8 (computed by this page's code at build time).
第 13 节的差分受步长两头夹击：截断误差要 h 小，舍入误差要 h 大。复数步长把这个两难拿掉：把一个坐标分量挪一个虚数步长 ih，整个 primal 求解用复数算，J 的虚部除以 h 就是导数：
The differences of Section 13 are squeezed from both sides: truncation wants a small h, round-off a large one. The complex step removes the dilemma: move one coordinate by an imaginary step ih, run the whole primal solve in complex arithmetic, and the imaginary part of J, divided by h, is the derivative:
原因是 Taylor 展开：
The reason is the Taylor expansion
实部与虚部各自分开：虚部的第一项就是 hdJ/dXm，下一项已是 h3。这里没有两个相近的数相减，也就没有相消，h 可以取 10−30：O(h2) 的截断误差远在舍入误差之下。
Real and imaginary parts separate: the imaginary part starts with hdJ/dXm, and its next term is already h3. No two nearly equal numbers are subtracted, so nothing cancels, and h can be 10−30: the O(h2) truncation error lies far below round-off.
把残差里的每一个实数换成复数：每个未知量、每个几何量都有实部与虚部两个数。挪动的是节点坐标，所以几何也变成复数：对偶面（边中点与单元中心）、面向量 n 与 |n|、边向量 ξ 与 p、k、最小二乘权重（其中 2×2 矩阵的逆）都要重算一遍。壁面行 ui−uw 里没有坐标，壁面值照旧直接写入。本页用到复数的加、减、乘、除，以及一次开方：|n|=n·n，取主值平方根。远场面是入流还是出流，沿用实数程序在基准网格上的判断，不随虚部变。两遍循环与 Jacobi 迭代（diag(Abc) 也由同一个复数循环收集）都不变，只是每个数多带一个虚部；判据要连虚部一起看：ImRbc/h 的均方根也要低于 10−13。从 u* 出发时实部一开始就满足判据，只看它，迭代在第 0 轮就会停下，漏掉状态的响应。伪代码在 图 14.1 的最上面。
Replace every real number in the residual by a complex one: every unknown and every geometric quantity carries a real and an imaginary part. What moves is a node coordinate, so the geometry becomes complex too: the dual faces (edge midpoints and cell centres), the face vector n and |n|, the edge vector ξ with p and k, and the least-squares weights (with the inverse of their 2×2 matrix) are all recomputed. The wall rows ui−uw contain no coordinate, and the wall values are written in as before. This page needs complex addition, subtraction, multiplication and division, and one square root: |n|=n·n, the principal root. Whether a far-field face is an inflow or an outflow face stays as the real code decided it on the base mesh; the imaginary part does not move it. The two loops and the Jacobi iteration (diag(Abc) collected by the same complex loop) stay as they are; every number just carries an imaginary part. The criterion must watch the imaginary part as well: the rms of ImRbc/h must fall below 10−13 too. Started at u*, the real part meets the criterion from the start, and a test on it alone would stop the iteration at sweep 0 and miss the state's response. The pseudocode is at the top of Figure 14.1.
复数迭代从哪里出发并非无关紧要。从收敛解 u* 出发（虚部为零），实部停在不动点上（到舍入误差为止）。把残差在 u* 处展开，Rbc(u*)=0：
Where the complex iteration starts matters. Started at the converged solution u* (imaginary part zero), the real part stays at the fixed point, to round-off. Expand the residual about u*, where Rbc(u*)=0:
于是 Jacobi 更新 u←u−diag(Abc)−1Rbc(u) 的虚部除以 h 是
so the imaginary part of the Jacobi update u←u−diag(Abc)−1Rbc(u), divided by h, is
这里 diag(Abc) 也是复数（它由同一个复数循环收集），所以更新里那个商的虚部完整写出来是
Here diag(Abc) is complex as well (the same complex loop collects it), so the imaginary part of the quotient in the update reads in full
从 u* 出发，ReRbc 停在舍入误差上，第二项也只有舍入误差的量级，剩下的正是第 10 节的前向 Jacobi 迭代：v0=0，同一个 diag(Abc)，右端项 f=∂Rbc/∂Xm（壁面行为 0）取自第 12 节的几何反向代码（逐行取出，是手写的导数，不是复数步长）。所以每一轮的 Imuk/h 就是前向迭代的 vk，不只是收敛后一致：48 个坐标分量、300 轮中的每一轮，最大差（除以 max|v|）是 2.2×10−15。实部与实数的 Jacobi 迭代相差不超过 1.1×10−16，但不是逐比特相同：实数程序用 Math.hypot 算 |n|，复数程序用平方根，两者在 80 个对偶面与半面中的 24 个上末位舍入不同。从初始流场出发就不同了：ReRbc≠0，第二项留下来，迭代把 diag(Abc) 也求了导，虚部是 primal 迭代值对坐标的导数：另一个序列，只在收敛时才到同一个极限。前几轮的差可达 max|v| 的 6.1 倍，而且每个分量都要把 primal 重解一遍：节点 1 的 x 从初始流场出发，虚部的残差要到第 127 轮才满足判据；从 u* 出发与前向迭代一样，是第 109 轮。
Started at u*, ReRbc sits at round-off, so the second term is down at round-off too, and what is left is the forward Jacobi iteration of Section 10: v0=0, the same diag(Abc), and the right-hand side f=∂Rbc/∂Xm (zero on the wall rows) taken from the geometric reverse code of Section 12 (row by row: a hand-written derivative, not a complex step). So Imuk/h is the forward iterate vk at every sweep, not only after convergence: over the 48 coordinates and each of the 300 sweeps the largest difference (over max|v|) is 2.2×10−15. The real part stays within 1.1×10−16 of the real Jacobi iterate, though not to the last bit: the real code takes |n| from Math.hypot and the complex code from a square root, and the two round the last bit differently on 24 of the 80 dual and half faces. Started at the initial field it is different: ReRbc≠0, the second term stays, the iteration differentiates diag(Abc) as well, and the imaginary part is the derivative of the primal iterate with respect to the coordinate: another sequence, which reaches the same limit only at convergence. In the first sweeps the two differ by up to 6.1 times max|v|, and every component re-solves the primal: for the x of node 1, started at the initial field, the imaginary residual meets the criterion at sweep 127; started at u*, at sweep 109, the forward iteration's count.
每个坐标分量一次复数求解，得到梯度的一个分量。48 个分量与伴随梯度之差的最大值（除以梯度的最大分量）是 3.0×10−16，第 13 节的中心差分（h=10−6）是 3.4×10−9：那里的差全是步长造成的。复数步长是精确的前向模式，所以它与伴随相符到舍入误差，是一个与步长无关的梯度检验。代价与前向模式一样，每个分量一次复数求解，而伴随一次线性求解就给出全部 48 个分量。在 图 14.1 里选一个坐标分量和起点，看两条迭代逐轮重合或分开。
One complex solve per coordinate gives one component of the gradient. Over the 48 components its largest difference from the adjoint gradient (over the gradient's largest component) is 3.0×10−16, against 3.4×10−9 for the central differences of Section 13 (h=10−6): that difference was all step size. The complex step is an exact forward mode, so it agrees with the adjoint to round-off, a gradient check whose result does not depend on the step size. It costs what forward mode costs, one complex solve per component, where one linear adjoint solve gives all 48. In Figure 14.1 pick a coordinate and a start, and watch the two iterations coincide sweep by sweep, or part.
伴随给出的是离散格式的精确导数，与格式本身是几阶精度无关；格式的精度阶决定的是离散解，以及由它算出的 J 与 dJ/dX，离连续问题有多远。这一节用构造解把它测出来：先选一个光滑函数当精确解，由它反推出源项与边界数据，再在一族逐次一致加密的网格上求解，看误差随网格尺寸按几次方减小。
The adjoint gives the exact derivative of the discrete scheme, whatever that scheme's order of accuracy; the order decides how far the discrete solution, and the J and dJ/dX computed from it, are from the continuous problem. This section measures it with a manufactured solution: pick a smooth function as the exact solution, derive the source term and the boundary data from it, solve on a family of uniformly refined meshes, and see with what power of the mesh size the error falls.
它处处介于 0.25 与 0.75 之间，与本页算例的解同号（u>0），所以对流速度 uβ·n 的符号分布与原问题相同。它的一阶与二阶导数（其中 0.325=14(0.92+0.72)）：
It lies between 0.25 and 0.75 everywhere, of the same sign as this page's solution (u>0), so the convective speed uβ·n has the same sign pattern as in the original problem. Its first and second derivatives (with 0.325=14(0.92+0.72)):
在方程右端加一个源项 s，让 ue 恰好是它的解：
Put a source term s on the right-hand side and make ue its exact solution:
这里用到 β 是常数：∇·(12u2β)=uβ·∇u。离散时每个控制体只取源项的一个值（单点积分），从残差里减去：
This uses β being constant: ∇·(12u2β)=uβ·∇u. The discrete source takes one value per control volume (a one-point quadrature) and is subtracted from the residual:
其中 |Vi| 是节点 i 的中位对偶控制体的面积。
where |Vi| is the area of node i's median-dual control volume.
壁面：uw=ue，照样强施加，壁面节点的行是 ui−ue(xi)。远场：入流半面的外侧状态取 ue，在半面中点 xf 求值；原问题里为零的远场扩散通量，换成精确解的法向通量：
Wall: uw=ue, imposed strongly as before, so the wall nodes' rows are ui−ue(xi). Far field: the outer state on an inflow half-face is ue at the half-face midpoint xf, and the far-field diffusive flux, zero in the original problem, becomes the exact solution's normal flux:
于是离散问题与连续问题只差在离散本身：把 ue 代入离散方程，剩下的残差全是截断误差。
The discrete problem now differs from the continuous one only by the discretisation itself: put ue into the discrete equations and what remains of the residual is truncation error, nothing else.
从 24 个单元的基础网格出发，每加密一层：每个三角形用三条边的中点分成 4 个三角形，每个四边形用四条边的中点与顶点平均点分成 4 个四边形，边界边在中点处一分为二。物面与远场都是直边八边形，所以每一层的计算域完全相同，网格尺寸逐层减半，hℓ=h0/2ℓ。六层网格有 24、96、384、1536、6144、24576 个单元，24、80、288、1088、4224、16640 个节点。拟合精度阶必须用这样一致细分出来的网格族；各自独立生成、只是目标尺寸不同的网格不构成网格族：网格的形状与误差常数 C 在各层之间都会变，拟合出来的“阶”没有保证。
Starting from the 24-cell base mesh, each refinement splits every triangle into 4 triangles through its edge midpoints and every quadrilateral into 4 quadrilaterals through its edge midpoints and its vertex average, and every boundary edge in two at its midpoint. The body and the far field are both straight-edged octagons, so the domain is exactly the same on every level and the mesh size halves from one level to the next, hℓ=h0/2ℓ. The six meshes have 24, 96, 384, 1536, 6144 and 24576 cells, and 24, 80, 288, 1088, 4224 and 16640 nodes. An order must be fitted on a family refined uniformly like this; meshes generated independently at different target sizes do not form a family: the mesh shape, and with it the error constant C, changes from one mesh to the next, so an "order" fitted over them carries no guarantee.
误差在节点上量，用对偶控制体面积加权的 L2 范数与最大范数：
The error is measured at the nodes, in the L2 norm weighted by the dual areas and in the max norm:
相邻两层的网格尺寸之比正好是 2，所以观测阶是误差之比以 2 为底的对数：
The mesh sizes of two neighbouring levels differ by exactly a factor of 2, so the observed order is the base-2 logarithm of the error ratio:
这里的 p 是观测阶，与第 3 节的两点系数 p=ξ·n/|ξ|2 无关。
Here p is the observed order, not the two-point coefficient p=ξ·n/|ξ|2 of Section 3.
每一层都用第 6 节的 Jacobi 迭代求解：先迭代到判据满足，比如在第 k 轮，再多跑 k 轮，这时已经落在舍入误差平台上。网格越密，Jacobi 越慢：最细一层只有扩散时要 61259 轮才满足判据，完整方程 39469 轮。从第 4 层起，本页代码本身跑得太慢，这几层用同一个 Jacobi 迭代的另一份写法来跑：把本页代码的每一步运算按原来的顺序搬到连续数组上，它在第 0–3 层与本页代码逐轮、逐位相同，只是快得多。离散解与求解器无关，误差表不受影响。本页格式在两个问题上的结果：
Every level is solved by the Jacobi iteration of Section 6: iterated until the criterion is met, at some sweep k, and then k sweeps more, by which time it is on its round-off floor. The finer the mesh, the slower Jacobi: on the finest level diffusion alone takes 61259 sweeps to meet the criterion, the full equation 39469. From level 4 on the page's own code is too slow, and those levels run a second transcription of the same Jacobi iteration: every operation of the page's code, in the same order, on flat arrays. On levels 0 to 3 it matches the page's code sweep for sweep and bit for bit, only much faster. The discrete solution does not depend on the solver, and the error table is not affected. This page's scheme on the two problems:
| 层 ℓ | level ℓ | 节点数 | nodes | 只有扩散（β = 0，ε = 0） | diffusion alone (β = 0, ε = 0) | 完整方程（β = (1, 0.5)，ε = 0.8） | the full equation (β = (1, 0.5), ε = 0.8) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ‖e‖2 | p | ‖e‖∞ | p | ‖e‖2 | p | ‖e‖∞ | p | ||||||||||||
| 0 | 24 | 1.83e−2 | — | 4.15e−2 | — | 3.89e−2 | — | 7.54e−2 | — | ||||||||||
| 1 | 80 | 4.63e−3 | 1.99 | 1.12e−2 | 1.90 | 2.69e−2 | 0.53 | 5.77e−2 | 0.39 | ||||||||||
| 2 | 288 | 1.09e−3 | 2.09 | 3.12e−3 | 1.84 | 1.65e−2 | 0.71 | 3.58e−2 | 0.69 | ||||||||||
| 3 | 1088 | 2.58e−4 | 2.08 | 7.97e−4 | 1.97 | 9.35e−3 | 0.82 | 2.05e−2 | 0.81 | ||||||||||
| 4 | 4224 | 6.23e−5 | 2.05 | 2.01e−4 | 1.99 | 5.02e−3 | 0.90 | 1.11e−2 | 0.88 | ||||||||||
| 5 | 16640 | 1.52e−5 | 2.03 | 5.00e−5 | 2.01 | 2.61e−3 | 0.94 | 5.80e−3 | 0.94 | ||||||||||
图 15.1 把两种格式在两个问题上的误差画在双对数坐标里：左边是解 u 的 L2 误差，右边是目标函数的误差（见下文“目标函数的精度阶”）。
Figure 15.1 plots the errors of both schemes on both problems on logarithmic axes: on the left the L2 error of u, on the right the error of the objective (see "The order of accuracy of the objective" below).
取 β=0、ε=0，方程只剩扩散与源项。本页 L2 范数下的阶是 1.99、2.09、2.08、2.05、2.03，最大范数在最细一层也到了 2.01；格心页面对应的 L2 阶是 1.94、2.02、2.01（格心格式只算到第 3 层：再加密，Jacobi 迭代在那里发散，见格心页面第 15 节）。二阶的前提是相容：最小二乘梯度与带修正的平均梯度通量对线性场都精确（第 3 节）。测得的阶说明，单点积分的源项与不规则的网格都没有把它拉低。
With β=0 and ε=0 only diffusion and the source remain. The L2 orders on this page are 1.99, 2.09, 2.08, 2.05 and 2.03, and the max-norm order reaches 2.01 on the finest level; the cell-centred page's L2 orders are 1.94, 2.02 and 2.01 (the cell-centred scheme only to level 3: refined further, the Jacobi iteration diverges there; see Section 15 of the cell-centred page). Behind second order is consistency: the least-squares gradients and the corrected average-gradient flux are exact on linear fields (Section 3). The measured orders show that neither the one-point source quadrature nor the irregular mesh pulls it down.
对流通量里的耗散项可以写成扩散两点部分的样子：
The dissipation in the convective flux can be written in the form of the two-point part of diffusion:
若 ξ 与 n 平行，ξ·n/|ξ|2=|n|/|ξ|，于是 νnum=ε|ξ|=O(h)：离散格式实际上是在解扩散系数为 ν+O(h) 的方程，误差是 O(h)，格式只有一阶。本页内部面上 νnum 的中位数逐层是 0.732、0.355、0.177、0.089、0.044、0.022，每加密一次大致减半；在基础网格上它比 ν=0.5 还大，到最细一层只剩 ν 的 4%。所以观测阶是从下面一步步靠近 1 的：L2 范数 0.53、0.71、0.82、0.90、0.94，最大范数 0.39、0.69、0.81、0.88、0.94（格心页面的 L2 阶是 0.59、0.71、0.82、0.90、0.95）。
If ξ is parallel to n, ξ·n/|ξ|2=|n|/|ξ| and so νnum=ε|ξ|=O(h): the discrete scheme is in effect solving the equation with diffusivity ν+O(h), the error is O(h), and the scheme is first order only. On this page the median of νnum over the interior faces is 0.732, 0.355, 0.177, 0.089, 0.044 and 0.022 level by level, roughly halving with every refinement; on the base mesh it exceeds ν=0.5, and on the finest it is 4% of ν. That is why the observed order approaches 1 from below, step by step: 0.53, 0.71, 0.82, 0.90 and 0.94 in L2, 0.39, 0.69, 0.81, 0.88 and 0.94 in the max norm (the cell-centred page's L2 orders are 0.59, 0.71, 0.82, 0.90 and 0.95).
上面量的是解 u 的误差；优化器用到的是目标函数 J，它的误差要单独量。构造解下 J 有精确值：精确解穿过物面 Γw 的总通量
What was measured above is the error of the solution u; what the optimiser uses is the objective J, and its error has to be measured on its own. Under the manufactured solution J has an exact value, the total flux of the exact solution through the body surface Γw,
（n̂ 为单位法向，指向物体内部）。物面是八段直线，每段用 24 点 Gauss–Legendre 求积，积分精确到舍入误差：只有扩散时 Je=−0.065629，完整方程 Je=−0.167841。离散的 Jh 就是本页第 7 节的目标函数，反力形式 J=−Σi∈WRi(u)，替换之前的壁面行里包括源项，在构造解的离散问题上求值：
(n̂ the unit normal, pointing into the body). The body surface is eight straight segments, and 24-point Gauss–Legendre quadrature on each gives the integral to round-off: Je=−0.065629 with diffusion alone, Je=−0.167841 for the full equation. The discrete Jh is this page's objective of Section 7, the reaction form J=−Σi∈WRi(u), the wall rows before the replacement including the source, evaluated on the discrete manufactured problem:
| 层 ℓ | level ℓ | 节点数 | nodes | 只有扩散（β = 0，ε = 0） | diffusion alone (β = 0, ε = 0) | 完整方程（β = (1, 0.5)，ε = 0.8） | the full equation (β = (1, 0.5), ε = 0.8) | ||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Jh − Je | p | Jh − Je | p | ||||||||
| 0 | 24 | −1.61e−2 | — | −7.15e−2 | — | ||||||
| 1 | 80 | −4.32e−3 | 1.90 | −3.46e−2 | 1.05 | ||||||
| 2 | 288 | −1.10e−3 | 1.98 | −1.72e−2 | 1.01 | ||||||
| 3 | 1088 | −2.75e−4 | 1.99 | −8.61e−3 | 1.00 | ||||||
| 4 | 4224 | −6.89e−5 | 2.00 | −4.31e−3 | 1.00 | ||||||
| 5 | 16640 | −1.72e−5 | 2.00 | −2.16e−3 | 1.00 | ||||||
J 的误差按格式的阶减小：只有扩散时观测阶是 1.90、1.98、1.99、2.00、2.00，完整方程是 1.05、1.01、1.00、1.00、1.00（格心页面：1.85、1.97、1.99；0.70、0.91、0.99、1.01、1.01）。完整方程里，u 的 L2 阶到最细一层是 0.94，还在从下面靠近 1；J 的阶在最细两层是 1.00、1.00，已经是 1。所以对这个目标函数，精度阶就是格式的阶：只有扩散时二阶，完整方程一阶。
The error of J falls at the scheme's order: the observed orders are 1.90, 1.98, 1.99, 2.00 and 2.00 with diffusion alone and 1.05, 1.01, 1.00, 1.00 and 1.00 for the full equation (on the cell-centred page: 1.85, 1.97 and 1.99; 0.70, 0.91, 0.99, 1.01 and 1.01). For the full equation the L2 order of u is 0.94 on the finest level, still approaching 1 from below, while the order of J on the two finest levels is 1.00 and 1.00, already 1. So for this objective the order of accuracy is the scheme's: second order with diffusion alone, first order for the full equation.
第 11–14 节验证的是：伴随与几何梯度对离散的 J 精确到舍入误差，这与格式几阶无关。离散的 J 本身随网格加密逼近连续值，速度由格式的阶决定：构造解下它的误差按格式的阶减小（上一小节）。本页自己的问题没有精确值，但两种格式的壁面通量在同一族网格上互相靠拢（本页第 7 节与格心页面第 7 节）：格心 2.344、2.331、2.249、2.184、2.145、2.124，格点 2.749、2.454、2.291、2.201、2.153、2.127。格点的相邻两层之差逐层按固定比例缩小，最后两个相邻之比是 1.85、1.91，接近一阶格式的 2。格心的相邻两层之差前几层忽大忽小，到最细几层才稳定：最后两个相邻之比是 1.66、1.83，正向一阶格式的 2 靠拢。用最细两层按一阶外推，格心指向 2.103，格点指向 2.102。优化器拿到的梯度与它所优化的离散目标函数完全一致，问题是自洽的；网格加密决定的是这个离散问题离连续问题有多近。要让它更近，得提高格式的阶（例如对流通量用重构），而不是改伴随。
What Sections 11 to 14 verified is that the adjoint and the geometric gradient are exact, to round-off, for the discrete J, whatever the scheme's order. The discrete J itself approaches the continuous value under refinement, at a pace set by the scheme's order: under the manufactured solution its error falls at the scheme's order (previous subsection). The page's own problem has no exact value, but on this same family the two schemes' wall fluxes close in on each other (Section 7 here and Section 7 of the cell-centred page): cell-centred 2.344, 2.331, 2.249, 2.184, 2.145 and 2.124, node-centred 2.749, 2.454, 2.291, 2.201, 2.153 and 2.127. The node-centred steps shrink by a steady factor, the last two ratios of successive steps being 1.85 and 1.91, close to the 2 of a first-order scheme. The cell-centred steps between levels are irregular at first and settle only on the finest levels: the last two ratios of successive steps are 1.66 and 1.83, moving towards the 2 of a first-order scheme. Extrapolated at first order from the two finest levels, the cell-centred values point to 2.103 and the node-centred ones to 2.102. The gradient the optimiser gets agrees exactly with the discrete objective it optimises, so the problem is self-consistent; refinement decides how close that discrete problem is to the continuous one. To bring it closer one has to raise the scheme's order (a reconstruction in the convective flux, say), not change the adjoint.
连续参考采用独立 P1 有限元，壁面伴随值为 1，远场按实现的总通量解释为 Robin 条件。格点强加行的乘子不直接代表壁面场值。加密与远场值梯度提供一致性证据，但不是连续形状梯度或任意网格收敛的证明。
The independent P1 finite-element reference has wall dual value one and a Robin far field derived from the implemented total flux. Strong node-row multipliers are not wall field samples. Refinement and the far-field-value gradient provide evidence, not a continuous shape-gradient or arbitrary-mesh convergence proof.
点积测试检查离散转置，差分和复数步长检查同一离散目标的梯度。独立连续参考回答另一问题：加密后是否趋向相同的连续灵敏度。误差不下降时，先核对边界、目标、乘子含义、正则性和参考误差，再讨论适定性。
Dot tests check the discrete transpose; finite differences and complex steps check the gradient of the same discrete objective. An independent continuous reference asks a different question: whether refinement approaches the same continuous sensitivity. A discrepancy calls for checking boundaries, objectives, multiplier interpretation, regularity and reference error before diagnosing well-posedness.
附录 A集中给出完整推导、参考实现、图与交互；可先继续阅读离散伴随主线。
Appendix A contains the full derivation, reference implementation, figures and interaction. You can continue the discrete-adjoint lesson first.
正文使用普通数组点积定义 。若另选带权内积，算子的伴随就会改变：状态权矩阵为 、残差权矩阵为 时，带权伴随是 。代码的普通转置不能在没有换元说明的情况下被写成该带权算子。
The main derivation defines using ordinary array dot products. A weighted inner product changes the adjoint: with state and residual weights and it is . The coded ordinary transpose cannot be relabelled as that weighted operator without a change of variables.
残差的归一化也会改变乘子的数值。若改为 ， 为可逆的行尺度矩阵（例如各控制体的体积），则在收敛解上 ，；目标总梯度保持不变。即使 随设计改变，由 产生的项也乘以 ，在精确收敛时为零。把控制体积分残差改成体积平均残差后，不能仍直接比较原来的乘子曲线。
Residual normalization also changes multiplier values. For with an invertible row-scaling matrix (for example the control-volume sizes), at convergence and , while the total objective gradient is unchanged. Even when depends on the design, the terms from multiply and vanish at exact convergence. Multipliers of integrated and volume-normalized residuals cannot be compared without rescaling.
连续附录从微分方程、积分目标和分部积分出发；正文从实际离散残差及其有限维链式法则出发。两者可以比较同一设计变量的目标总导数与加密极限，但还须匹配符号、变量、积分权重和边界约束。被替换的边界行的乘子尤其不能直接解释为连续场采样值。
The continuous appendix starts from differential equations, integral objectives and integration by parts; the main text starts from actual discrete residuals and finite-dimensional chain rules. Compare total derivatives of the same objective and design variable and their refinement limits, after matching signs, variables, weights and boundary constraints. Multipliers of replaced boundary rows especially are not direct samples of a continuous field.
| 本页（格点） | this page (node-centred) | 格心页面 | cell-centred page | ||
|---|---|---|---|---|---|
| 壁面条件 | wall condition | 强施加：壁面节点的行换成 ui−uw | strong: the wall nodes' rows become ui−uw | 弱施加：uw 进入壁面面的通量 | weak: uw enters the wall faces' flux |
| 梯度模板 | gradient stencil | 边邻居 | edge neighbours | 面邻居，外加壁面面的中点（取 uw） | face neighbours, plus the wall faces' midpoints (carrying uw) |
| 目标函数 | objective | 反力：−Σi∈WRi(u) | reaction: −Σi∈WRi(u) | 壁面面通量之和 | the sum of the wall-face fluxes |
| 伴随里的壁面 | the wall in the adjoint | 循环前清零、循环后加回 | zero before the loop, add back after | 没有特殊处理：壁面面就是一个普通的边界面 | nothing special: a wall face is an ordinary boundary face |
| Jacobi 的谱半径 ρ(I − diag(A)−1A)，格点用 Abc | Jacobi spectral radius ρ(I − diag(A)−1A), with Abc on the node-centred page | 0.787 | 0.915 | ||
| 所算最细一层的 L2 阶（只有扩散 / 完整方程） | L2 order on the finest level computed (diffusion alone / full equation) | 2.03 / 0.94 | 2.01 (ℓ = 3) / 0.95 | ||
| 所算最细一层目标函数误差的阶（只有扩散 / 完整方程） | order of the objective's error on the finest level computed (diffusion alone / full equation) | 2.00 / 1.00 | 1.99 (ℓ = 3) / 1.01 |
本项目以离散伴随为主线：从可逐句检查的小模型，逐步走向工业级非结构网格 RANS 的完整离散伴随实现。连续伴随在简单算例中提供独立验证与对比，不是学习后续章节的前提。
This project follows discrete adjoints, from small models whose operations can be inspected statement by statement toward a complete industrial unstructured-grid RANS discrete adjoint. Continuous adjoints provide independent checks and comparisons for simple cases; they are not a prerequisite for later chapters.
这里的“走向”是项目目标，并不表示当前已实现 RANS。可以为选定的连续 RANS 模型推导形式伴随；但不能据此假定它严格对应工业求解器中的湍流闭合、壁面处理、限幅、边界与网格算法。最终需要验证的是实际离散残差及目标的导数。
RANS is the project destination, not a capability already implemented here. A formal continuous adjoint can be derived for a specified continuous RANS model; this does not establish strict equivalence to an industrial solver’s turbulence closure, wall treatment, limiters, boundary and mesh algorithms. The derivatives to verify are those of the implemented discrete residual and objective.
两条独立计算路径，在统一的目标、边界与误差范数下比较
Two independent paths compared with matched objectives, boundaries and norms
连续 PDE → 线性化与分部积分 → 连续边界条件 → 独立参考场 λ
Continuous PDE → linearize and integrate by parts → dual boundary conditions → independent λ
离散残差 Rh → 精确 Jacobian → 转置求解 → 守恒行上的 −ψh
Discrete residual Rh → exact Jacobian → transpose solve → −ψh on conservation rows
将定常 PDE 写为 ，目标泛函记为 ，采用拉格朗日泛函
页面求解的是 。这里的连续乘子采用相反符号：在守恒方程上有 。这些残差已经是控制体积分后的守恒平衡，不需要再除以控制体的体积或面积。被代数边界约束替换的行对应约束乘子，其数值不是连续场的点值。
三类检查分别回答不同的问题：
有限组网格上误差下降是数值证据，不是收敛定理。即使积分范数下的误差下降，边界层或特征跳跃仍可能阻止最大范数下的收敛。
出现差异提供了诊断机会，但差异本身不能证明问题不适定。应先明确状态空间与边界迹空间、允许的数据、目标泛函以及范数。Hadamard 意义下的适定性要求解存在、唯一且连续依赖于数据。正确推导的连续伴随可以是适定的，而离散边界闭合却不一致；一个边界数据与目标泛函不匹配的适定方程，仍然是错误的伴随。约束乘子、有限正则性和参考解误差也会造成表观差异。
把原始边界条件的齐次线性化约束纳入算子 的定义域 。如果该算子是适当空间之间的有界同构，则其 Banach 伴随也是相应对偶空间之间的有界同构。这给出一种充分的局部分析框架，并不是写出形式微分表达式就能得到的结论。对于无界 PDE 算子，必须先明确其定义域及图范数，或者采用适当的弱形式。对于非线性原始问题，仅有存在唯一性还不能保证解映射可微；需要在所选解支上检查局部光滑性及线性化可逆性。
对于 ，原始算子的零模态给出必要的兼容性条件：
若 含有非零元素，则伴随一旦有解就不唯一。仅有与原始零空间正交的条件，并不能构成一般的无穷维存在性定理：还需要考虑值域闭性和稳定性。对齐次伴随边界数据，所需估计具有以下形式：
如果边界数据非齐次，还应加上相应的边界迹范数。每个网格尺度下逆算子有界也不够：在正确加权的离散范数下，稳定界必须与网格尺度无关。不考虑积分残差和控制体权重的原始欧氏条件数，本身不能作为连续问题的诊断结论。
取 ，固定原始入口值，目标变分为 。对原始线性化 使用 Green 恒等式，得到
对流扩散。 正扩散系数提供椭圆主部，但这不是任意低阶项与 Robin 边界系数下的普遍可逆性证明。仍需检查齐次零空间，以及所用能量估计或最大值原理的假设。独立计算支持当前参数组，没有证明对所有参数都成立的定理。把目标所要求的壁面伴随值 1 改成 0，可能得到另一个适定问题，但会给出错误的目标灵敏度。在远场同时规定状态值和法向导数，一般会过约束。把被替换方程行的乘子当成壁面场值，则是另一类代数含义错误。
因此，教学顺序应当是：先明确连续原始问题及所选解支；推导伴随，包括其边界定义域；分析存在性、唯一性和稳定性；随后才检查对该问题的一致、稳定逼近。反例应明确保留其失败的前提，不能悄悄替换模型后宣布验证成功。
此时有
对状态的变分 ，完整的 Green 恒等式是
域内伴随中没有 ，因为原始线性化采用守恒形式。
壁面采用固定 Dirichlet 数据 。目标是沿流体域外法向流入物体的总通量：
固定几何时，壁面上 ，但 可以自由变化。因此
将其与 的边界项合并，得到 。于是
对这个通量目标，令壁面伴随值为零是错误的。
令 ，其中 是单位外法向。将积分后的数值面通量除以面长，代码实现的远场总通量为
相容的连续解释只有一个条件：。特别地，在入流处，它是非线性 Robin 条件：
代码中的“扩散面通量为零”与外部对流状态，不能被理解为同时施加连续条件 和 。在出流处， 才简化为 。
边界线性化给出 ，因此远场的自由变分导出
注意 仍然保留在这个边值问题中。作为内部人工扩散，它的影响随网格加密消失，但单位长度的边界通量仍然含有这一系数。
对于入流数据，驻值条件给出一个可独立核对的标量灵敏度：
格点页采用 ，并将壁面方程替换为 。把原始 Jacobian 按自由节点与壁面节点分块。伴随的自由部分满足
因此应当比较的连续场是 ，并取 ，因为
原始的 是边界约束乘子。把它画成连续壁面迹，会人为制造表观不一致。
为完整起见，令 表示光滑区域变形，并保持物理参数 为常数。入流／出流划分不变，没有面跨越 。记 。在驻值状态，记对流通量 ，拉格朗日泛函化为
由于 ，壁面目标与残差的壁面边界项相消。保持拉回参考域后的状态与乘子不变，对该表达式求导，得到
其中入流处 ，出流处 。如果随变形输运的边界数据不是常数，还需要加入其物质导数。这里给出了该公式的推导，但新增实验实际验证的是 ，不是完整形状导数；原页面的坐标导数检查仍然是离散检查。
P1 有限元参考在同一个多边形区域的三角剖分上求解非线性原始弱式，然后直接组装推导得到的连续伴随弱式
并在壁面规定 。该过程不调用有限体积转置。三角形与边上的积分规则用于计算多项式被积函数。参考网格分别有 1,088、4,224 和 16,640 个顶点，用于检查参考分辨率。最终参考仍然是数值近似。离线研究中的有限体积状态方程和伴随方程使用原页面的数值核心，通过稀疏直接 Newton／线性求解将代数残差压低。这改变的是验证时的求解器，不是有限体积方程，也不是页面的 Jacobi 演示。
Write the steady PDE as , the functional as , and use
The pages solve . Our continuous multiplier has the opposite sign: on conservation equations. Those residuals are integrated balances, so there is no additional division by control-volume measure. Rows replaced by algebraic boundary constraints carry constraint multipliers; their values are not point samples of the continuous field.
Three tests answer different questions:
A decreasing error on finitely many meshes is numerical evidence, not a theorem of convergence. Boundary layers or characteristic jumps can prevent maximum-norm convergence even when an integral norm decreases.
A discrepancy is a diagnostic opportunity, not by itself proof of ill-posedness. Specify the state and trace spaces, admissible data, functional, and norms first. Hadamard well-posedness requires existence, uniqueness and continuous dependence. A correctly derived continuous dual can be well posed while a discrete boundary closure is inconsistent; a well-posed equation with the wrong objective boundary data is simply the wrong dual. Constraint multipliers, limited regularity and finite reference error are other sources of apparent disagreement.
Include the homogeneous linearized primal boundary constraints in the domain of the operator . If this is a bounded isomorphism between suitable spaces, its Banach adjoint is also a bounded isomorphism between the dual spaces. This is a sufficient local framework, not a conclusion obtained merely by writing a formal differential expression. An unbounded PDE operator must first be given a domain/graph norm or a suitable weak realization. For a nonlinear primal, existence and uniqueness alone do not guarantee a differentiable solution map; the local smoothness and invertibility hypotheses must be checked on the chosen branch.
For , a primal null mode implies the necessary compatibility condition
A nonzero element of makes a solvable dual nonunique. Orthogonality to the primal kernel alone is not a general infinite-dimensional existence theorem: closed range and stability still matter. The required estimate has the form
for homogeneous dual boundary data; inhomogeneous boundary data contribute their appropriate trace norm. A bounded inverse at every mesh size is not enough either: the stability bound must be uniform in the properly weighted discrete norms. Raw Euclidean condition numbers without accounting for integrated residuals and control-volume weights are not themselves a continuum diagnosis.
For , fixed primal inflow and objective variation , the Green identity for the primal linearization gives
Advection–diffusion. Positive diffusion provides an elliptic principal part, not a universal invertibility proof for arbitrary lower-order and Robin coefficients. Homogeneous kernels and the assumptions of applicable energy or maximum-principle estimates still need attention. The independent solves support this parameter set; they do not prove a theorem for all parameters. Setting the wall dual to zero instead of the objective-required one can define a different well-posed problem, but yields the wrong functional sensitivity. Imposing both state and normal derivative at the far field is generally an overconstraint. Treating a replaced-row multiplier as a wall field sample is yet another, algebraic interpretation error.
The teaching sequence is therefore: identify the continuous primal and chosen branch; derive the dual including its boundary domain; examine existence, uniqueness and stability; only then test consistent, stable approximation of that problem. A counterexample should retain its failed premise visibly, rather than silently substituting another model and declaring success.
Here
For a state variation , the complete Green identity is
There is no in the volume dual: the primal linearization is in conservative form.
The wall has fixed Dirichlet data . Its objective is the total outward fluid flux into the body,
At fixed geometry, on the wall but is free. Hence
Combining this with the boundary term of gives . Therefore
Zero wall adjoint data would be wrong for this flux functional.
Let for the unit outward normal. Dividing the integrated numerical face flux by face length, the implemented far-field total flux is
The consistent continuum interpretation is the single condition . In particular, on inflow this is a nonlinear Robin condition,
The code's “zero diffusive face flux” and external convective state must not be read as simultaneously imposing continuum and . On outflow, reduces to .
Linearization gives , so the free far-field variation yields
Notice that remains in this boundary-value problem. It vanishes with refinement as an interior artificial diffusion, but the boundary flux per unit length retains it.
For the incoming datum, stationarity gives an independently checkable scalar sensitivity:
The node page has and replaced rows . Partition the raw Jacobian into free and wall blocks. The free part of the adjoint satisfies
Thus the continuous field to compare is and , since
The raw is a boundary-constraint multiplier. Plotting it as the wall trace would manufacture an apparent inconsistency.
For completeness, let be a smooth domain deformation, with constant physical . The inflow/outflow partition is held fixed and no face crosses . Write . At a stationary state, with the convective flux , the Lagrangian reduces to
The wall objective cancels the wall residual boundary term because . Differentiating this expression while holding the pulled-back state and multiplier fixed gives
where on inflow and on outflow. Nonconstant transported boundary data would add their material derivatives. This formula is derived here but the added experiment validates , not this full shape derivative; the earlier pages' coordinate checks remain discrete checks.
A P1 finite-element reference solves the nonlinear primal weak form on a triangulation of the same polygon, then directly assembles the derived continuous dual weak form
with on the wall. It does not call the finite-volume transpose. Triangle and edge quadrature integrate the polynomial integrands. References on 1,088, 4,224 and 16,640 vertices check reference resolution. The final reference is still a numerical approximation. FV state and dual equations use the original page cores, solved to tight algebraic residuals with sparse direct Newton/linear solves for this offline study. This changes the verification solver, not the finite-volume equations or the page's Jacobi demonstrations.
三角剖分与形函数梯度 → Newton 求原始弱式 → 直接组装连续伴随弱式 → 施加 λ壁=1 → 求解并积分远场梯度
Triangulate and form basis gradients → Newton solve of primal weak form → assemble continuous dual weak form → impose wall λ=1 → solve and integrate far-field gradient
下方代码可以完整复制到本地文件运行。它独立计算连续参考，不调用本页有限体积 Jacobian 的转置。二维代码内嵌两次细化的多边形网格；这里输出参考场，后面的多网格交互保留已计算的加密结果。
Copy the complete code below into a local file and run it. It computes the continuous reference independently of the page’s finite-volume Jacobian transpose. The 2-D scripts include the polygonal mesh refined twice; they output reference fields, while the refinement viewer below retains the precomputed multi-grid comparisons.
python3 -m pip install numpy scipy
python3 reference_ad.pyimport numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
from scipy.integrate import solve_bvp,simpson
from scipy.optimize import brentq
from scipy.linalg import null_space
def triangulate(mesh):
X=np.array(mesh['X']);T=[]
for c in mesh['C2N']:
c=np.array(c)-1;T.extend([c[:3]] if len(c)==3 else [c[[0,1,2]],c[[0,2,3]]])
return X,np.array(T)
def fem(mesh,uinf=1.):
"""P1 Galerkin primal and directly assembled continuous dual on a fixed polygon.
Three-point triangle / three-point edge quadrature; no FV transpose used.
"""
X,T=triangulate(mesh);nt=len(T);nn=len(X);beta=np.array([1.,.5]);nu=.5;eps=.8
P=X[T];a=P[:,1]-P[:,0];b=P[:,2]-P[:,0];det=a[:,0]*b[:,1]-a[:,1]*b[:,0];V=det/2
assert min(V)>0
G=np.stack([np.stack([P[:,1,1]-P[:,2,1],P[:,2,0]-P[:,1,0]],axis=1),np.stack([P[:,2,1]-P[:,0,1],P[:,0,0]-P[:,2,0]],axis=1),np.stack([P[:,0,1]-P[:,1,1],P[:,1,0]-P[:,0,0]],axis=1)],axis=1)/det[:,None,None]
bg=G@beta;K=nu*V[:,None,None]*np.einsum('tik,tjk->tij',G,G)
bary=np.array([[2/3,1/6,1/6],[1/6,2/3,1/6],[1/6,1/6,2/3]])
wall=set();edges=[]
# Directed CCW boundary edges: normal points out of the computational domain.
for c in mesh['C2N']:
for a,b in zip(c,c[1:]+c[:1]):
key=f'{min(a,b)}-{max(a,b)}'
if key not in mesh['bnd']:continue
if mesh['bnd'][key]=='wall':wall.update([a-1,b-1])
else:edges.append([a-1,b-1])
E=np.array(edges);d=X[E[:,1]]-X[E[:,0]];length=np.linalg.norm(d,axis=1);bn=(d[:,1]-.5*d[:,0])/length
zz,ww=np.polynomial.legendre.leggauss(3);edgephi=np.stack([(1-zz)/2,(1+zz)/2],axis=1);ww=ww/2
rows=np.broadcast_to(T[:,:,None],(nt,3,3)).ravel();cols=np.broadcast_to(T[:,None,:],(nt,3,3)).ravel()
erows=np.broadcast_to(E[:,:,None],(len(E),2,2)).ravel();ecols=np.broadcast_to(E[:,None,:],(len(E),2,2)).ravel()
fixed=np.array(sorted(wall));free=np.setdiff1d(np.arange(nn),fixed)
def assemble(u,dual=False):
ut=u[T];uq=ut@bary.T
local=K.copy();r=np.einsum('tij,tj->ti',K,ut)
for k,phi in enumerate(bary):
r-=V[:,None]/3*.5*uq[:,k,None]**2*bg
if dual:local-=V[:,None,None]/3*uq[:,k,None,None]*phi[None,:,None]*bg[:,None,:]
else:local-=V[:,None,None]/3*uq[:,k,None,None]*bg[:,:,None]*phi[None,None,:]
be=np.zeros((len(E),2,2));br=np.zeros((len(E),2))
for phi,w in zip(edgephi,ww):
q=u[E]@phi
h=np.where(bn<0,.25*bn*(q*q+uinf*uinf)+eps*(q-uinf),.5*bn*q*q)
hp=np.where(bn<0,.5*bn*q+eps,bn*q)
br+=length[:,None]*w*h[:,None]*phi
be+=length[:,None,None]*w*hp[:,None,None]*np.outer(phi,phi)
A=coo_matrix((np.r_[local.ravel(),be.ravel()],(np.r_[rows,erows],np.r_[cols,ecols])),shape=(nn,nn)).tocsr()
R=np.bincount(T.ravel(),weights=r.ravel(),minlength=nn)+np.bincount(E.ravel(),weights=br.ravel(),minlength=nn)
return R,A
u=np.ones(nn);u[fixed]=0
for it in range(30):
R,A=assemble(u)
if max(abs(R[free]))<2e-12:break
du=spsolve(A[free][:,free],-R[free]);u[free]+=du
else:raise RuntimeError('FEM primal not converged')
_,D=assemble(u,True);lam=np.ones(nn);lam[free]=spsolve(D[free][:,free],-D[free][:,fixed]@np.ones(len(fixed)))
dualres=float(max(abs((D@lam)[free])))
# Incoming datum sensitivity and conservative reaction flux.
grad=0.
for phi,w in zip(edgephi,ww):grad-=np.sum(length*w*(lam[E]@phi)*np.where(bn<0,.5*bn*uinf-eps,0))
# Same weak conservation law: wall reaction is minus assembled residual there.
J=-np.sum(R[fixed])
return dict(X=X,T=T,u=u,lam=lam,J=float(J),gradient=float(grad),primal_res=float(max(abs(R[free]))),dual_res=dualres)
import json
mesh=json.loads('{"X":[[0.804,-0.04],[0.6125777772109133,0.4725777772109134],[0.10000000000000003,0.51],[-0.165239682094289,0.22523968209428905],[-0.29600000000000004,-0.03999999999999995],[-0.25857777721091335,-0.3985777772109132],[0.09999999999999991,-0.5900000000000001],[0.519239682094289,-0.45923968209428906],[1.2972000000000001,-0.022000000000000002],[0.9414940753805006,0.8644940753805004],[0.05500000000000008,1.1355],[-0.695458123066357,0.7284581230663572],[-1.0178,-0.02199999999999987],[-0.7467940753805007,-0.8237940753805003],[0.05499999999999979,-1.1795],[0.8901581230663569,-0.8571581230663573],[1.9,0],[1.3435028842544403,1.34350288425444],[1.1634144591899854e-16,1.9],[-1.34350288425444,1.3435028842544403],[-1.9,2.326828918379971e-16],[-1.3435028842544405,-1.34350288425444],[-3.490243377569956e-16,-1.9],[1.3435028842544399,-1.3435028842544405],[1.0506000000000002,-0.031],[1.1193470376902503,0.4212470376902502],[0.777035926295707,0.6685359262957069],[0.7082888886054567,0.2162888886054567],[0.9138179631478536,0.31876796314785344],[0.4982470376902503,0.9999970376902502],[0.07750000000000005,0.82275],[0.3562888886054567,0.4912888886054567],[0.42726796314785354,0.7456429631478534],[-0.3202290615331785,0.9319790615331786],[-0.430348902580323,0.4768489025803231],[-0.03261984104714448,0.3676198410471445],[-0.17642445129016146,0.6497994512901616],[-0.8566290615331785,0.3532290615331786],[-0.6569,-0.03099999999999991],[-0.2306198410471445,0.09261984104714455],[-0.5436244512901616,0.22292445129016158],[-0.8822970376902504,-0.4228970376902501],[-0.5026859262957071,-0.6111859262957068],[-0.2772888886054567,-0.2192888886054566],[-0.5797929631478536,-0.32109296314785335],[-0.34589703769025043,-1.0016470376902502],[0.07749999999999985,-0.88475],[-0.07928888860545671,-0.4942888886054566],[-0.21259296314785361,-0.7479679631478533],[0.47257906153317836,-1.0183290615331786],[0.704698902580323,-0.6581989025803232],[0.3096198410471444,-0.5246198410471445],[0.3910994512901613,-0.7714744512901617],[1.0936790615331784,-0.43957906153317866],[0.6616198410471446,-0.24961984104714452],[0.8776494512901616,-0.3445994512901616],[1.5986,-0.011000000000000001],[1.62175144212722,0.67175144212722],[1.3203514421272202,0.66075144212722],[1.1424984798174704,1.1039984798174702],[0.6717514421272203,1.6217514421272199],[0.47074703769025034,1.38224703769025],[0.027500000000000097,1.51775],[-0.6717514421272199,1.62175144212722],[-0.64425144212722,1.2395014421272201],[-1.0194805036603984,1.0359805036603986],[-1.6217514421272199,0.6717514421272203],[-1.2977290615331785,0.3642290615331787],[-1.4588999999999999,-0.010999999999999819],[-1.6217514421272203,-0.6717514421272199],[-1.1806514421272203,-0.6827514421272199],[-1.0451484798174706,-1.0836484798174701],[-0.6717514421272205,-1.6217514421272199],[-0.3733970376902505,-1.3618970376902502],[0.027499999999999723,-1.53975],[0.6717514421272197,-1.6217514421272203],[0.6992514421272198,-1.2615014421272202],[1.1168305036603985,-1.100330503660399],[1.6217514421272199,-0.6717514421272203],[1.3950790615331785,-0.42857906153317865],[0.9273000000000001,-0.035500000000000004],[0.9822089815739269,0.14388398157392673],[0.8110534258766551,0.26752842587665504],[0.7561444443027283,0.08814444430272834],[0.8691767129383277,0.11601421293832753],[1.1739000000000002,-0.026500000000000003],[1.2082735188451252,0.19962351884512508],[1.016582500419052,0.37000750041905184],[1.0952412502095261,0.1717537502095259],[1.0304205565353755,0.6428705565353753],[0.8592650008381038,0.7665150008381036],[0.8454269447217803,0.4936519447217802],[0.9379237506285779,0.5682612506285778],[0.6948068517533101,0.5705568517533102],[0.6604333329081851,0.344433332908185],[0.7529301388149826,0.4190426388149826],[0.6021519447217802,0.7070894447217801],[0.3917784258766551,0.618465925876655],[0.484433332908185,0.4819333329081851],[0.5432926388149826,0.5945113888149827],[0.7198705565353755,0.9322455565353753],[0.46275750041905195,0.8728200004190518],[0.6610112506285778,0.8196675006285776],[0.2766235188451252,1.067748518845125],[0.06625000000000006,0.979125],[0.2523839815739268,0.7841964815739266],[0.264503750209526,0.9259725002095259],[0.08875000000000005,0.6663749999999999],[0.22814444430272837,0.5006444443027284],[0.24026421293832756,0.6424204629383274],[-0.0494622256450807,0.7362747256450808],[-0.10452214616865296,0.5087096461686531],[0.03369007947642778,0.43880992052357226],[-0.00788607308432646,0.5875423230843264],[-0.1326145307665892,1.0337395307665893],[-0.24832675641166996,0.7908892564116701],[-0.09103837820583496,0.885007128205835],[-0.5078435922997677,0.8302185922997678],[-0.56290351282334,0.6026535128233401],[-0.30338667693524224,0.5633241769352424],[-0.405615134617505,0.696771384617505],[-0.297794292337306,0.35104429233730605],[-0.09892976157071673,0.29642976157071677],[-0.20115821925297947,0.4298769692529796],[-0.4869866769352423,0.34988667693524234],[-0.38712214616865304,0.15777214616865307],[-0.19792976157071673,0.1589297615707168],[-0.34245821925297953,0.2544082192529796],[-0.7760435922997677,0.5408435922997679],[-0.70012675641167,0.2880767564116701],[-0.631515134617505,0.4453651346175051],[-0.9372145307665893,0.16561453076658939],[-0.83735,-0.02649999999999989],[-0.6002622256450808,0.09596222564508083],[-0.768738378205835,0.13078837820583508],[-0.47645000000000004,-0.035499999999999934],[-0.26330992052357227,0.026309920523572297],[-0.43178607308432654,0.061136073084326566],[-0.6183464815739268,-0.17604648157392663],[-0.42854092587665515,-0.27019092587665494],[-0.2866444443027284,-0.12964444430272826],[-0.4524954629383276,-0.15284546293832746],[-0.9500485188451252,-0.22244851884512498],[-0.731045000419052,-0.37199500041905176],[-0.784197500209526,-0.1992475002095258],[-0.8145455565353755,-0.6233455565353752],[-0.6247400008381039,-0.7174900008381035],[-0.5412394447217803,-0.46613944472178004],[-0.6778925006285779,-0.5447425006285777],[-0.3806318517533102,-0.50488185175331],[-0.267933332908185,-0.3089333329081849],[-0.40458638881498266,-0.3875363888149825],[-0.35763944472178033,-0.6795769447217801],[-0.14594092587665516,-0.621128425876655],[-0.16893333290818502,-0.4464333329081849],[-0.2632863888149827,-0.5630051388149825],[-0.5463455565353755,-0.9127205565353753],[-0.27924500041905204,-0.8748075004190518],[-0.45199250062857793,-0.7961487506285776],[-0.14544851884512533,-1.090573518845125],[0.06624999999999982,-1.032125],[-0.06754648157392688,-0.8163589815739267],[-0.1064975002095261,-0.9534662502095258],[0.08874999999999988,-0.7373750000000001],[0.010355555697271597,-0.5421444443027283],[-0.028595462938327643,-0.6792517129383275],[0.2342997256450806,-0.8281122256450808],[0.3503596461686529,-0.6480471461686531],[0.20480992052357216,-0.5573099205235723],[0.2195548230843264,-0.6927110730843267],[0.26378953076658906,-1.0989145307665893],[0.43183925641166987,-0.8949017564116701],[0.24904462820583484,-0.963513378205835],[0.6813685922997677,-0.9377435922997679],[0.7974285128233399,-0.7576785128233403],[0.5478991769352421,-0.7148366769352424],[0.6146338846175049,-0.8262901346175051],[0.6119692923373059,-0.5587192923373061],[0.41442976157071665,-0.4919297615707168],[0.4811644692529794,-0.6033832192529797],[0.7911741769352423,-0.5013991769352424],[0.7696346461686531,-0.2971096461686531],[0.5904297615707168,-0.3544297615707168],[0.6908019692529794,-0.4279144692529796],[0.9919185922997676,-0.648368592299768],[0.98566425641167,-0.39208925641167014],[0.8915463846175049,-0.5748838846175053],[1.1954395307665893,-0.23078953076658934],[0.9641247256450809,-0.1877997256450808],[1.079782128205835,-0.20929462820583508],[0.7328099205235723,-0.14480992052357225],[0.8484673230843267,-0.16630482308432656],[1.4479000000000002,-0.0165],[1.45947572106361,0.32487572106361],[1.3087757210636102,0.31937572106361],[1.7492999999999999,-0.0055000000000000005],[1.7608757210636101,0.33587572106361],[1.61017572106361,0.33037572106361],[1.4710514421272203,0.6662514421272201],[1.48262716319083,1.00762716319083],[1.3319271631908303,1.00212716319083],[1.2198492399087353,0.5409992399087351],[1.2430006820359554,1.2237506820359552],[1.2314249609723453,0.8823749609723451],[1.1309227587538604,0.7626227587538602],[1.0419962775989855,0.9842462775989853],[0.8066227587538604,1.2431227587538602],[0.7061205565353754,1.123370556535375],[1.0076271631908302,1.48262716319083],[0.9071249609723453,1.362874960972345],[0.5712492399087353,1.501999239908735],[0.3358757210636102,1.76087572106361],[0.23537351884512522,1.641123518845125],[0.4844970376902503,1.19112203769025],[0.013750000000000108,1.708875],[0.2491235188451252,1.449998518845125],[0.2628735188451252,1.258873518845125],[0.041250000000000085,1.326625],[-0.30837572106360994,1.3786257210636101],[-0.29462572106360996,1.1875007210636102],[-0.3358757210636099,1.7608757210636101],[-0.3221257210636099,1.5697507210636101],[-0.65800144212722,1.4306264421272201],[-1.00762716319083,1.48262716319083],[-0.9938771631908301,1.29150216319083],[-0.48224025183019925,1.0857402518301993],[-1.1814916939574194,1.1897416939574195],[-0.8318659728938091,1.1377409728938095],[-0.6698547825967884,0.9839797825967886],[-0.8574693133633777,0.8822193133633779],[-1.1586047825967885,0.7001047825967887],[-0.9965935922997677,0.5463435922997679],[-1.48262716319083,1.0076271631908302],[-1.3206159728938092,0.8538659728938094],[-1.4597402518301992,0.5179902518301995],[-1.76087572106361,0.33587572106361024],[-1.5988645307665892,0.18211453076658946],[-1.0771790615331784,0.35872906153317863],[-1.6794499999999999,-0.005499999999999793],[-1.3783145307665892,0.17661453076658942],[-1.1577645307665891,0.1711145307665894],[-1.23835,-0.016499999999999845],[-1.31977572106361,-0.34687572106360987],[-1.0992257210636103,-0.3523757210636099],[-1.7608757210636101,-0.33587572106360986],[-1.54032572106361,-0.34137572106360986],[-1.4012014421272203,-0.67725144212722],[-1.4826271631908305,-1.00762716319083],[-1.2620771631908303,-1.01312716319083],[-1.0314742399087353,-0.5528242399087351],[-1.1943256820359556,-1.2135756820359551],[-1.1128999609723453,-0.8831999609723451],[-0.9637227587538605,-0.7532727587538601],[-0.8959712775989856,-0.9537212775989852],[-0.7092727587538605,-1.2227727587538602],[-0.5600955565353756,-1.0928455565353752],[-1.0076271631908305,-1.48262716319083],[-0.8584499609723455,-1.3526999609723451],[-0.5225742399087355,-1.491824239908735],[-0.3358757210636104,-1.76087572106361],[-0.18669851884512542,-1.6309485188451252],[-0.35964703769025047,-1.18177203769025],[0.013749999999999686,-1.719875],[-0.17294851884512538,-1.450823518845125],[-0.15919851884512534,-1.270698518845125],[0.04124999999999976,-1.3596249999999999],[0.36337572106360977,-1.40062572106361],[0.3771257210636098,-1.22050072106361],[0.3358757210636097,-1.7608757210636101],[0.34962572106360973,-1.5807507210636103],[0.6855014421272198,-1.4416264421272202],[1.0076271631908298,-1.4826271631908305],[1.02137716319083,-1.3025021631908302],[0.5859152518301991,-1.1399152518301994],[1.2301666939574192,-1.2219166939574198],[0.9080409728938091,-1.1809159728938097],[0.7947047825967885,-1.0593297825967887],[1.0034943133633778,-0.9787443133633782],[1.2559547825967885,-0.7644547825967888],[1.1426185922997676,-0.642868592299768],[1.4826271631908299,-1.0076271631908305],[1.3692909728938092,-0.8860409728938097],[1.5084152518301992,-0.5501652518301995],[1.76087572106361,-0.33587572106361013],[1.6475395307665892,-0.21428953076658933],[1.2443790615331785,-0.43407906153317866],[1.4968395307665894,-0.21978953076658933],[1.3461395307665893,-0.22528953076658934]],"C2N":[[1,81,85,84],[81,25,82,85],[85,82,29,83],[84,85,83,28],[25,86,89,82],[86,9,87,89],[89,87,26,88],[82,89,88,29],[29,88,93,92],[88,26,90,93],[93,90,10,91],[92,93,91,27],[28,83,96,95],[83,29,92,96],[96,92,27,94],[95,96,94,2],[2,94,100,99],[94,27,97,100],[100,97,33,98],[99,100,98,32],[27,91,103,97],[91,10,101,103],[103,101,30,102],[97,103,102,33],[33,102,107,106],[102,30,104,107],[107,104,11,105],[106,107,105,31],[32,98,110,109],[98,33,106,110],[110,106,31,108],[109,110,108,3],[3,108,114,113],[108,31,111,114],[114,111,37,112],[113,114,112,36],[31,105,117,111],[105,11,115,117],[117,115,34,116],[111,117,116,37],[37,116,121,120],[116,34,118,121],[121,118,12,119],[120,121,119,35],[36,112,124,123],[112,37,120,124],[124,120,35,122],[123,124,122,4],[4,122,128,127],[122,35,125,128],[128,125,41,126],[127,128,126,40],[35,119,131,125],[119,12,129,131],[131,129,38,130],[125,131,130,41],[41,130,135,134],[130,38,132,135],[135,132,13,133],[134,135,133,39],[40,126,138,137],[126,41,134,138],[138,134,39,136],[137,138,136,5],[5,136,142,141],[136,39,139,142],[142,139,45,140],[141,142,140,44],[39,133,145,139],[133,13,143,145],[145,143,42,144],[139,145,144,45],[45,144,149,148],[144,42,146,149],[149,146,14,147],[148,149,147,43],[44,140,152,151],[140,45,148,152],[152,148,43,150],[151,152,150,6],[6,150,156,155],[150,43,153,156],[156,153,49,154],[155,156,154,48],[43,147,159,153],[147,14,157,159],[159,157,46,158],[153,159,158,49],[49,158,163,162],[158,46,160,163],[163,160,15,161],[162,163,161,47],[48,154,166,165],[154,49,162,166],[166,162,47,164],[165,166,164,7],[7,164,170,169],[164,47,167,170],[170,167,53,168],[169,170,168,52],[47,161,173,167],[161,15,171,173],[173,171,50,172],[167,173,172,53],[53,172,177,176],[172,50,174,177],[177,174,16,175],[176,177,175,51],[52,168,180,179],[168,53,176,180],[180,176,51,178],[179,180,178,8],[8,178,184,183],[178,51,181,184],[184,181,56,182],[183,184,182,55],[51,175,187,181],[175,16,185,187],[187,185,54,186],[181,187,186,56],[56,186,190,189],[186,54,188,190],[190,188,9,86],[189,190,86,25],[55,182,192,191],[182,56,189,192],[192,189,25,81],[191,192,81,1],[9,193,195],[193,57,194],[195,194,59],[193,194,195],[57,196,198],[196,17,197],[198,197,58],[196,197,198],[59,199,201],[199,58,200],[201,200,18],[199,200,201],[57,198,194],[198,58,199],[194,199,59],[198,199,194],[9,195,87],[195,59,202],[87,202,26],[195,202,87],[59,201,204],[201,18,203],[204,203,60],[201,203,204],[26,205,90],[205,60,206],[90,206,10],[205,206,90],[59,204,202],[204,60,205],[202,205,26],[204,205,202],[10,206,208],[206,60,207],[208,207,62],[206,207,208],[60,203,210],[203,18,209],[210,209,61],[203,209,210],[62,211,213],[211,61,212],[213,212,19],[211,212,213],[60,210,207],[210,61,211],[207,211,62],[210,211,207],[10,208,101],[208,62,214],[101,214,30],[208,214,101],[62,213,216],[213,19,215],[216,215,63],[213,215,216],[30,217,104],[217,63,218],[104,218,11],[217,218,104],[62,216,214],[216,63,217],[214,217,30],[216,217,214],[11,218,220],[218,63,219],[220,219,65],[218,219,220],[63,215,222],[215,19,221],[222,221,64],[215,221,222],[65,223,225],[223,64,224],[225,224,20],[223,224,225],[63,222,219],[222,64,223],[219,223,65],[222,223,219],[11,220,115],[220,65,226],[115,226,34],[220,226,115],[65,225,228],[225,20,227],[228,227,66],[225,227,228],[34,229,118],[229,66,230],[118,230,12],[229,230,118],[65,228,226],[228,66,229],[226,229,34],[228,229,226],[12,230,232],[230,66,231],[232,231,68],[230,231,232],[66,227,234],[227,20,233],[234,233,67],[227,233,234],[68,235,237],[235,67,236],[237,236,21],[235,236,237],[66,234,231],[234,67,235],[231,235,68],[234,235,231],[12,232,129],[232,68,238],[129,238,38],[232,238,129],[68,237,240],[237,21,239],[240,239,69],[237,239,240],[38,241,132],[241,69,242],[132,242,13],[241,242,132],[68,240,238],[240,69,241],[238,241,38],[240,241,238],[13,242,244],[242,69,243],[244,243,71],[242,243,244],[69,239,246],[239,21,245],[246,245,70],[239,245,246],[71,247,249],[247,70,248],[249,248,22],[247,248,249],[69,246,243],[246,70,247],[243,247,71],[246,247,243],[13,244,143],[244,71,250],[143,250,42],[244,250,143],[71,249,252],[249,22,251],[252,251,72],[249,251,252],[42,253,146],[253,72,254],[146,254,14],[253,254,146],[71,252,250],[252,72,253],[250,253,42],[252,253,250],[14,254,256],[254,72,255],[256,255,74],[254,255,256],[72,251,258],[251,22,257],[258,257,73],[251,257,258],[74,259,261],[259,73,260],[261,260,23],[259,260,261],[72,258,255],[258,73,259],[255,259,74],[258,259,255],[14,256,157],[256,74,262],[157,262,46],[256,262,157],[74,261,264],[261,23,263],[264,263,75],[261,263,264],[46,265,160],[265,75,266],[160,266,15],[265,266,160],[74,264,262],[264,75,265],[262,265,46],[264,265,262],[15,266,268],[266,75,267],[268,267,77],[266,267,268],[75,263,270],[263,23,269],[270,269,76],[263,269,270],[77,271,273],[271,76,272],[273,272,24],[271,272,273],[75,270,267],[270,76,271],[267,271,77],[270,271,267],[15,268,171],[268,77,274],[171,274,50],[268,274,171],[77,273,276],[273,24,275],[276,275,78],[273,275,276],[50,277,174],[277,78,278],[174,278,16],[277,278,174],[77,276,274],[276,78,277],[274,277,50],[276,277,274],[16,278,280],[278,78,279],[280,279,80],[278,279,280],[78,275,282],[275,24,281],[282,281,79],[275,281,282],[80,283,285],[283,79,284],[285,284,17],[283,284,285],[78,282,279],[282,79,283],[279,283,80],[282,283,279],[16,280,185],[280,80,286],[185,286,54],[280,286,185],[80,285,287],[285,17,196],[287,196,57],[285,196,287],[54,288,188],[288,57,193],[188,193,9],[288,193,188],[80,287,286],[287,57,288],[286,288,54],[287,288,286]],"bnd":{"1-84":"wall","28-84":"wall","2-95":"wall","28-95":"wall","17-197":"far","58-197":"far","18-200":"far","58-200":"far","2-99":"wall","32-99":"wall","3-109":"wall","32-109":"wall","18-209":"far","61-209":"far","19-212":"far","61-212":"far","3-113":"wall","36-113":"wall","4-123":"wall","36-123":"wall","19-221":"far","64-221":"far","20-224":"far","64-224":"far","4-127":"wall","40-127":"wall","5-137":"wall","40-137":"wall","20-233":"far","67-233":"far","21-236":"far","67-236":"far","5-141":"wall","44-141":"wall","6-151":"wall","44-151":"wall","21-245":"far","70-245":"far","22-248":"far","70-248":"far","6-155":"wall","48-155":"wall","7-165":"wall","48-165":"wall","22-257":"far","73-257":"far","23-260":"far","73-260":"far","7-169":"wall","52-169":"wall","8-179":"wall","52-179":"wall","23-269":"far","76-269":"far","24-272":"far","76-272":"far","1-191":"wall","55-191":"wall","8-183":"wall","55-183":"wall","17-284":"far","79-284":"far","24-281":"far","79-281":"far"}}')
if __name__ == '__main__':
r=fem(mesh);h=1e-4
fd=(fem(mesh,1+h)['J']-fem(mesh,1-h)['J'])/(2*h)
np.savetxt('continuous_ad.csv',np.column_stack([r['X'],r['u'],r['lam']]),delimiter=',',header='x,y,u,lambda')
print('primal residual',r['primal_res'],'dual residual',r['dual_res'])
print('continuous gradient',r['gradient'],'finite difference',fd)
点积测试查转置，差分与复数步长查同一个离散目标的梯度；它们都说明不了离散伴随是否逼近连续伴随。本页的连续伴随问题在附录 A.1 第 3 小节推导：壁面上 λ=1，远场是由实现的总通量得到的 Robin 条件，以及 dJ/du∞ 的连续公式。这里只给出一个独立参考与本页有限体积离散在加密下的比较。比较时离散乘子按 λh=−ψh 换号（附录 A.1 第 1 小节）；守恒残差已经对控制体积分，不再除以体积。
A dot test checks the transpose, and finite differences and the complex step check the gradient of the same discrete objective; none of them says whether the discrete adjoint approaches the continuous one. This page's continuous adjoint problem is derived in Appendix A.1, part 3: λ=1 on the wall, a Robin condition at the far field that comes from the implemented total flux, and the continuous formula for dJ/du∞. This part only compares an independent reference with the page's finite-volume discretisation under refinement. The discrete multipliers are compared after the sign change λh=−ψh (Appendix A.1, part 1); the conservation residuals are already integrated over the control volumes and are not divided by the volume again.
本页的壁面行是约束行：比较时壁面上的场取重建值 λW=1（附录 A.1 第 3 小节），原始强加行上的 −ψW 是约束乘子，不能拿来与连续壁面值比较；下面的表、下图及误差都用这个重建场，代数残差检查仍保留所有未知量。
The wall rows of this page are constraint rows: for the comparison the wall field takes the reconstructed value λW=1 (Appendix A.1, part 3), and the raw −ψW on the strongly imposed rows are constraint multipliers, not samples of the continuous wall value. The table below, the viewer and the errors all use the reconstructed field; algebraic residual checks keep every unknown.
在同一个多边形区域上用 P1 有限元独立求解原始 PDE，再直接组装附录 A.1 第 3 小节推导的连续伴随弱式求解；不调用有限体积转置。参考网格 1,088、4,224、16,640 个顶点。最细连续梯度为 3.292e+00；4,224 顶点的伴随梯度 3.293e+00 与重解原始 PDE 的中心差分 3.293e+00 一致。参考仍有离散误差，不能称作精确解析解。
On the same polygon, an independent P1 finite-element primal solve is followed by direct assembly of the continuous dual weak form derived in Appendix A.1, part 3, without calling the FV transpose. Reference meshes have 1,088, 4,224 and 16,640 vertices. The finest reference gradient is 3.292e+00; the 4,224-vertex dual gradient 3.293e+00 agrees with the independently recomputed PDE central difference 3.293e+00. The reference still has discretisation error; it is not an exact analytic solution.
| 层Level | DOF | L2 | 参考差Ref. gap | dJ/du∞ |
|---|---|---|---|---|
| 0 | 24 | 1.038e-01 | 4.931e-04 | 3.690e+00 |
| 1 | 80 | 6.559e-02 | 4.781e-04 | 3.513e+00 |
| 2 | 288 | 4.091e-02 | 4.793e-04 | 3.421e+00 |
| 3 | 1088 | 2.377e-02 | 4.819e-04 | 3.366e+00 |
L2 按控制体面积加权。参考差是两级最细独立参考场在同一取样点上的 L2 差，比这里的有限体积误差小得多。加密给出积分范数下一致性的证据，不是最大范数或任意网格上的收敛证明；边界和角点还可能限制收敛率。
L2 is weighted by control-volume area. The reference gap compares the two finest independent references at the same sample points and is much smaller than the FV discrepancy. Refinement provides evidence in this integral norm, not a maximum-norm or arbitrary-mesh theorem; boundaries and corners may limit rates.
表中有限体积各层的解，解的是本页原来的残差与精确 Jacobian；这张表离线用稀疏 Newton／直接线性求解计算，离散解与求解器无关。这里数值验证的只是远场值的导数 dJ/du∞；附录 A.1 第 3 小节写出的完整形状变分，没有在这里做数值验证。
The finite-volume levels of the table solve this page's own residual with its exact Jacobian; the table was computed offline with sparse Newton / direct linear solves, and the discrete solution does not depend on the solver. What is verified numerically here is only the derivative with respect to the far-field value, dJ/du∞; the full shape variation written in Appendix A.1, part 3, is not verified numerically here.
dJ/du∞ 这一列。表中的值从上方单调靠近参考值 3.292。表只到第 3 层；用本页自己的 Jacobi 迭代（第 10 节的切线求解，每层先到判据、再跑同样多轮）把同一族网格算到第 5 层，得到 3.690、3.513、3.421、3.366、3.332、3.313，前 4 层与表中一致，之后继续单调下降。
The dJ/du∞ column. The table's values approach the reference 3.292 monotonically from above. The table stops at level 3; the page's own Jacobi iteration (the tangent solve of Section 10, every level run to the criterion and as many sweeps again) carries the same family to level 5: 3.690, 3.513, 3.421, 3.366, 3.332 and 3.313, the first 4 as in the table, still falling monotonically after that.
结果对不上时应该先查什么（存在性、唯一性、稳定性，边界条件与目标是否匹配，乘子的含义），在附录 A.1 第 2 小节（选读）。
What to check first when results disagree (existence, uniqueness and stability, whether the boundary conditions and the objective match, what the multipliers mean) is in Appendix A.1, part 2 (optional reading).
下图展示已离线计算的连续参考与离散场，可切换网格；它不是在浏览器里重新求解连续 PDE。图 A.1 给出场与误差。完整的符号、边界项、形状变分和证据范围见连续伴随推导。
The viewer displays offline-computed continuous references and discrete fields at selectable resolutions; it does not rerun a continuous PDE solve in the browser. Figure A.1 shows fields and errors. See the continuous-adjoint derivation for signs, boundary terms, shape variations and evidence limits.