本页的方程是对流方程的格心页面上的方程加一项常系数扩散 −ν∇2u（ν=0.5），算例是同一张环形（O 型）混合网格。有了扩散，内圈物面才能施加 Dirichlet 壁面条件 u=uw=0。格心格式的未知量在单元里，物面上没有未知量，所以壁面值只能弱施加：它进入壁面面上的通量，也作为一个点进入壁面单元的最小二乘梯度。扩散通量要用单元梯度，残差因此成了两遍循环：先算梯度，再算通量。目标函数是流进物面的总通量。对照阅读：同一方程的格点页面，那里的壁面是强施加的。
The equation is that of the cell-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. A cell-centred scheme keeps its unknowns in the cells and has none on the body, so the wall value can only be imposed weakly: it enters the flux on the wall faces, and it is one more point in the wall cells' least-squares gradients. The diffusive flux needs cell gradients, so the residual becomes two loops: gradients first, then fluxes. The objective is the total flux into the body. Read alongside: the node-centred page for the same equation, where the wall is imposed strongly.
把 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:
第一项是对流页面的非线性守恒律，第二项是常系数扩散。求解域 Ω 是物体与远场之间的环形区域。源项并没有消失：每个单元一个源项强度，作为设计变量留在离散方程里，第 11 节对它求导，第 15 节用它构造精确解。
The first term is the nonlinear conservation law of the advection page, the second constant-coefficient diffusion. The domain Ω is the annulus between the body and the far field. The source has not gone away: one source strength per cell 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* 在物面附近最小，在迎着来流的远场附近最大；下面几小节逐一说明这张图。
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 smallest near 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 把基准解上的两种通量逐点画出来。按每一圈的总和看，贴壁的一圈里扩散通量占上风，最外一圈里对流通量占上风；贴壁一圈的扩散箭头都偏向物体。
Figure 1.2 draws the two fluxes of the baseline solution point by point. Summed ring by ring, diffusion wins in the ring against the wall and convection in the outermost ring; in the ring against the wall every diffusive arrow leans towards the body.
物体是一个不规则的八边形，远场是半径 1.9 的圆的内接八边形，两者之间就是求解域，与对流页面相同。数据是 β=(1,0.5)、ν=0.5、来流值 u∞=1 与壁面值 uw=0。基准解 u* 在 24 个单元上介于 0.179 与 0.879 之间：最小的在背着来流一侧的贴壁单元 c1，最大的在迎着来流的远场单元 c17（图 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. Over the 24 cells the baseline solution u* lies between 0.179 and 0.879: the smallest in the wall cell c1 on the side facing away from the flow, the largest in the far-field cell c17 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，就是 8 个壁面面上的通量之和。在基准解上 J=2.343758：从远场的 4 个入流面流入 2.569307，从 4 个出流面流出 0.225549，两者之差全部流进物体。壁面上 u=0，对流通量 12u2β·n 为零，所以流进物体的全是扩散通量。
The objective of Section 7 is the total flux into the body, J, the sum of the fluxes through the 8 wall faces. At the baseline solution J=2.343758: 2.569307 flows in through the 4 inflow faces of the far field and 0.225549 out through the 4 outflow 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∞,sc,uw,f,X) 有四部分：来流值 u∞；每个单元 c 里的源项强度 sc（基准值为 0）；每个壁面面 f 上的壁面值 uw,f（基准值为 0）；以及 24 个网格节点的坐标 X（网格见第 2 节）。合起来是 1 + 24 + 8 + 48 = 81 个分量。按伴随优化文献的通常写法，设计变量记作 D。第 10、11 节求 J 对 u∞ 的导数，第 11 节求它对 sc 与 uw,f 的导数，第 12 节求它对 X 的导数。
The design input D=(u∞,sc,uw,f,X) has four parts: the free-stream value u∞; a source strength sc in each cell c (0 in the baseline); a wall value uw,f on each wall face f (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 sc and uw,f, 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 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 in the cells; diffusion also needs two new geometric quantities and the gradient stencil of each cell.
网格与对流页面第 2 节完全相同：24 个节点，贴壁一层是 8 个四边形，外层是 16 个三角形，一共 24 个单元。未知量放在单元里，控制体就是单元本身。每条网格边就是一个面：32 个内部面，8 个壁面面，8 个远场面，一共 48 个面。
The mesh is exactly the one of Section 2 of the advection page: 24 nodes, with a layer of 8 quadrilaterals against the wall and 16 triangles outside it, 24 cells in all. The unknowns sit in the cells, and each control volume is the cell itself. Every mesh edge is a face: 32 interior faces, 8 wall faces and 8 far-field faces, 48 faces in all.
物面上没有未知量：它只是 8 个壁面面，各自属于一个贴壁的四边形单元（图 2.1）。壁面值 uw 因此不能像格点格式那样直接写进某个未知量，只能通过这 8 个面进入方程（第 4 节）。
There are no unknowns on the body: it is just the 8 wall faces, each belonging to one quadrilateral against the wall (Figure 2.1). So the wall value uw cannot be written into an unknown, as the node-centred scheme does; it can only enter the equations through those 8 faces (Section 4).
每个单元的未知量代表单元里的一个点：本页取它各顶点的平均点 xc。三角形的顶点平均点就是形心；四边形上两者一般不同，本网格 8 个贴壁四边形上，两点的距离是 |Ωc| 的 3.5% 到 8.6%（|Ωc| 为单元面积）。用顶点平均点有两个原因。第一，它是顶点坐标的线性函数，第 12 节的几何导数里，它得到的反向权重正好平均分给单元的各个顶点（三角形每个 13，四边形每个 14）。第二，取哪一点不影响格式对线性场的精确性：最小二乘梯度和修正后的扩散通量（第 3 节）只要求用到的值是线性场在所用点上的值。一致加密时四边形越来越接近平行四边形，上面那个比例的最大值在第 15 节的各层网格上依次是 8.6%、6.6%、4.4%、2.6%、1.4%、0.7%：它随网格尺寸一起趋于零，最后几层每层减半。
The unknown of a cell stands for one point in the cell, and this page takes the average of its vertices, xc. On a triangle the vertex average is the centroid; on a quadrilateral the two differ in general, and on the 8 quadrilaterals against the wall of this mesh they are 3.5% to 8.6% of |Ωc| apart (|Ωc| the cell's area). There are two reasons for the vertex average. First, it is linear in the vertex coordinates: in the geometric derivative of Section 12 the reverse weight it receives is shared equally among the cell's vertices (13 each on a triangle, 14 on a quadrilateral). Second, the choice of point does not affect the scheme's exactness on linear fields: the least-squares gradient and the corrected diffusive flux (Section 3) only need the values to be the linear field's values at the points used. Under uniform refinement the quadrilaterals approach parallelograms, and the largest value of the ratio above is 8.6%, 6.6%, 4.4%, 2.6%, 1.4% and 0.7% on the levels of Section 15: it goes to zero with the mesh size, halving per level on the last levels.
扩散项要两样新的几何量。第一样是面两侧未知量之间的向量 ξ：内部面上从左单元的 xL 指向右单元的 xR，壁面面上从壁面单元的 xL 指向面中点（壁面值 uw 就放在那里）。第二样是每个单元的梯度模板：它的面邻居的顶点平均点，再加上它自己壁面面的中点（如果有的话）。远场面不带数据（扩散通量为零），不进入模板。于是壁面单元有 4 个模板点（3 个面邻居加 1 个壁面中点），内层三角形有 3 个，外层三角形只有 2 个：两个点在二维里正好确定一个梯度。
Diffusion needs two new geometric quantities. The first is the vector ξ between the unknowns on either side of a face: on an interior face from the left cell's xL to the right cell's xR, on a wall face from the wall cell's xL to the face midpoint (where the wall value uw sits). The second is each cell's gradient stencil: the vertex averages of its face neighbours, plus the midpoint of its own wall face if it has one. Far-field faces carry no data (their diffusive flux is zero) and stay out of the stencils. So a wall cell has 4 stencil points (3 face neighbours and 1 wall midpoint), an inner triangle 3, and an outer triangle only 2: in two dimensions two points fix a gradient exactly.
在 图 2.2 里点选任意单元，可以看到它的面和梯度模板（虚线箭头，琥珀色的指向壁面中点）；再点选一个面，可以看到它的 n、ξ，以及下一节要用的分解 n=pξ+k。
In Figure 2.2, pick any cell to see its faces and its gradient stencil (the dashed arrows, amber for a wall midpoint); then pick a 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 cell and writes the convective and the diffusive flux on every face, with the least-squares gradients the diffusive flux needs.
本节写出每个面上的通量：先把方程化成各个面上通量之和（有限体积离散），再写对流通量（与对流页面相同）和扩散通量；扩散通量要面上的法向导数，因此还要最小二乘梯度。
This section writes the flux on every face: first the equation becomes a sum of 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 节的方程在单元 c 上积分，用散度定理化成边界积分：
Integrate the equation of Section 1 over cell c and turn the volume integral into a boundary integral with the divergence theorem:
求和取遍单元 c 的所有面，包括它的壁面面与远场面。每个面上的通量只算一次，以相反符号分给左右两个单元（左取 +、右取 −），格式是守恒的。ha 是对流通量（下标 a 取自 advection），hv 是扩散通量。这里的 n 是面的法向（从 L 指向 R，边界面上指向域外）乘以面长，|n| 就是面长，所以 vn、ha、hv 都已含面长；上面积分里的 n 则是单位法向。有源项时，平衡式再减去源项在单元上的积分：Rc←Rc−sc；基准问题中 sc=0，它作为设计变量留在方程里（第 11 节）。
The sum runs over all faces of cell c, its wall and far-field faces included. Each face's flux is computed once and handed to its left and right cells with opposite signs (+ left, − right), so the scheme is conservative. ha is the convective flux (the subscript is from advection) and hv the diffusive one. Here n is the face normal (from L to R, out of the domain on a boundary face) scaled by the face length, so |n| is the face length and vn, ha and hv already carry it; in the integral above n is the unit normal. With a source, the balance also subtracts the source integrated over the cell, Rc←Rc−sc; sc=0 in the baseline, and it stays in the equations as a design variable (Section 11).
中心平均加常数耗散，远场上按冻结的入流、出流取右状态，推导在对流页面第 3、4 节。壁面面是例外，那里状态已知，见第 4 节。
A central average plus constant dissipation, with the right state on the far field taken by the frozen inflow/outflow switch, derived in Sections 3 and 4 of the advection page. The wall faces are the exception, since the state there is known; see Section 4.
扩散通量 hv=−ν∇u·n 要面上的法向导数，但面两侧只有两个单元值。只有当 n 与两个顶点平均点的连线 ξ 平行时，两点差才就是法向导数；一般网格上并不平行（图 3.1）。于是把 n 分成沿 ξ 的部分和垂直于 ξ 的部分：
The diffusive flux hv=−ν∇u·n needs the normal derivative on the face, but only two cell values sit either side of it. The two-point difference is the normal derivative only when n is parallel to the vector ξ joining the two vertex averages, which on a general mesh it is not (Figure 3.1). So n is split into a part along ξ and a part perpendicular to it:
沿 ξ 的部分用两点差，垂直的部分用两侧单元梯度的平均。两点部分让左右两个单元直接耦合，梯度部分补上 n 偏离 ξ 的那一块；两部分对线性场都精确。
The part along ξ takes the two-point difference, the perpendicular part the average of the two cell gradients. The two-point part couples the left and right cells directly; the gradient part makes up for n leaning away from ξ. Both are exact on linear fields.
单元 c 的梯度 Gc 取为让线性模型 uc+G·ξs 最贴近各模板点上的值 us 的那个 G；模板点是面邻居的顶点平均点，或壁面面的中点（那里的值是 uw）：
The gradient Gc of cell c is the G for which the linear model uc+G·ξs comes closest to the values us at the stencil points, which are the vertex averages of the face neighbours, or the midpoint of a wall face (where the value is uw):
权重只依赖几何，线性场被精确复现（ΣswcsξsT=I），本页不加距离权。壁面中点让壁面值进入了壁面单元的梯度，这是弱施加的第二个入口。梯度对 u 是线性的，对 uw 也是。对状态求导时 uw 是常数，没有它的份（第 8 节）；把壁值本身当作设计变量求导，见第 11 节。图 3.3 可以逐个单元查看。
The weights depend on the geometry only, linear fields come out exactly (ΣswcsξsT=I), and there is no distance weighting. The wall midpoints bring the wall value into the wall cells' gradients: the second door through which the wall is imposed weakly. The gradient is linear in u, and in uw too. Differentiating with respect to the state, uw is a constant with no share (Section 8); for the derivative with respect to the wall value itself, see Section 11. Figure 3.3 shows it cell by cell.
上面的公式是一个最小二乘问题的解。下面从问题出发，一步一步推到公式，再推到第 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.
第一步，写出要最小化的量。在单元 c 的每个模板点 s 上，线性模型 uc+G·ξs 与那里的值 us 差 us−uc−G·ξs；最小二乘取使这些差的平方和最小的 G：
Step 1, the quantity to minimise. At each stencil point s of cell c the linear model uc+G·ξs misses the value us there by us−uc−G·ξs; least squares takes the G that makes the sum of their squares smallest:
第二步，令它对 G 的导数为零。平方和对 G 求导得 −2Σsξs(us−uc−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−uc−G·ξs); with ξs(G·ξs)=(ξsξsT)G this gives the 2×2 normal equations:
第三步，解出 Gc。左边的 2×2 矩阵 Mc（正文的记号）只依赖几何。由 Cauchy–Schwarz 不等式，它的行列式不小于零，只有当所有 ξs 都平行时才等于零（所以外圈三角形只有 2 个模板点也够）。两边乘它的逆，就是正文的公式：
Step 3, solve for Gc. The 2×2 matrix Mc 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 (so the 2 points of an outer triangle are enough). 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−uc=ξs·∇u=ξsT∇u，代入第三步：
Step 4, linear fields come out exactly. Let u be a linear field with a constant gradient ∇u; then us−uc=ξs·∇u=ξsT∇u, and step 3 gives
这就是 ΣswcsξsT=Mc−1Mc=I 的意思。壁面中点也在内：只要那里的 uw 也取自同一个线性场。
That is what ΣswcsξsT=Mc−1Mc=I says. The wall midpoints are included, provided the uw there comes from the same linear field.
第五步，对状态线性化。Gc 对 u 是线性的，偏导数逐项读出（m 不在模板里时 wcm=0）：
Step 5, the linearisation in the state. Gc is linear in u, and its partial derivatives can be read off term by term (wcm=0 when m is not in the stencil):
壁面中点上的 uw 是数据，不是未知量，没有自己的列；但它的权重照样计入 −Σswcs。这就是第 6 节 ∂Gc/∂us=wcs、∂Gc/∂uc=−Σswcs 的来历，也是代码 gradRows() 存下的量。
The uw at a wall midpoint is data, not an unknown, and has no column of its own; its weight still counts in −Σswcs. This is where the ∂Gc/∂us=wcs and ∂Gc/∂uc=−Σswcs of Section 6 come from, and what the code stores in gradRows().
第六步，转置：伴随怎样把第一遍倒过来。（记号：加点的 ẋ 是前向（tangent）模式的导数，加横线的 x̄ 是反向（adjoint）模式累加的伴随，第 8 节的“记号”一段有完整说明。）梯度的伴随 Ḡc 与梯度的切线配成下面的对偶恒等式；按状态的切线 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 Ḡc 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 cell itself gives one up:
一个内部面在两个模板里各出现一次（R 在 L 的模板里，L 在 R 的模板里）。按面组织这个循环，就是 residual_b 第 10 至 14 行的 γL=wLR·ḠL、γR=wRL·ḠR（代码把 Ḡ 存在数组 G_b 里，ū 存在 y 里）。壁面面只在 L 的模板里出现，而它的 s 是壁面中点，没有 ūs 可写：那一份丢弃（第 12 行只对内部面执行）。
An interior face appears in two stencils (R in that of L, L in that of R). Organised by faces, this loop is lines 10 to 14 of residual_b, γL=wLR·ḠL and γR=wRL·ḠR (the code keeps Ḡ in the array G_b and ū in y). A wall face appears only in the stencil of L, and its s is the wall midpoint, with no ūs to write: that share is dropped (line 12 runs for interior faces only).
核对：本页 24 个单元上，按第三步的闭式算出的权重与代码的 wcs 最多差 5.8×10−16；ΣswcsξsT−I 的最大分量是 5.6×10−16；线性场 u=0.3+0.7x−0.4y 的最小二乘梯度与 (0.7,−0.4) 最多差 1.1×10−15；gradRows() 给出的 ∂G/∂u 与对 grads() 的中心差分最多差 5.0×10−14（本页代码在构建时计算）。
Check: over the 24 cells of this page, the weights from the closed form of step 3 differ from the code's wcs by at most 5.8×10−16; the largest entry of ΣswcsξsT−I is 5.6×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 1.1×10−15; and the ∂G/∂u of gradRows() differs from central differences of grads() by at most 5.0×10−14 (computed by this page's code at build time).
图 3.2 画出三种单元的模板向量 ξs 与权重 wcs：权重一般不与 ξs 同向，正是 2×2 矩阵的逆把它们调成对线性场精确的组合。
Figure 3.2 draws the stencil vectors ξs and the weights wcs of three kinds of cell: 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.
图 3.4 把一个内部面或壁面面上的量全部列出来，状态取收敛解 u*，并给出 ∇u·n 的三种估计：只用平均梯度、只用两点差，以及修正后的两者之和。
Figure 3.4 lists everything on one interior or wall face at the converged solution u*, with three estimates of ∇u·n: the averaged gradient alone, the two-point difference alone, and their corrected sum.
第一步，p 由“k 与 ξ 垂直”唯一确定：
Step 1, p is fixed by requiring k to be perpendicular to ξ:
第二步，对线性场精确。设 u 是梯度为常向量 ∇u 的线性场。ξ 连接 uL、uR 所在的两点，所以 uR−uL=ξ·∇u；由上一段第四步，两侧的最小二乘梯度都等于 ∇u，所以平均梯度 Gf=∇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 Gf=∇u. Hence
两点部分和梯度部分各管 n 的一段，加起来正好是 n。壁面面上 Gf=GL、uR=uw，只要 uw 取自同一个线性场，结论不变。
The two-point part and the gradient part each take one piece of n, and together they take all of it. On a wall face Gf=GL and uR=uw, and the conclusion stands as long as uw comes from the same linear field.
第三步，对四个输入线性化：两侧的状态 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:
壁面面上 Gf=GL，uR=uw 是数据，对流部分 12uw2vn 与 uL 无关（cL=0）：
On a wall face Gf=GL, uR=uw is data, and the convective part 12uw2vn does not depend on uL (cL=0):
第四步，这些偏导数就是 residual_d 第 11 行 ḣ=cLu̇L+cRu̇R−ν[p(u̇R−u̇L)+Ġf·k] 里的系数，也是 residual_b 反过来乘的系数。再经第一遍链式求导（上一段第五步），就得到第 6 节的公式（[m=L] 在 m=L 时为 1，否则为 0）：
Step 4, these partial derivatives are the coefficients on line 11 of residual_d, ḣ=cLu̇L+cRu̇R−ν[p(u̇R−u̇L)+Ġf·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):
核对：在本页 40 个内部面与壁面面上，对 faceFlux() 取中心差分（步长 10−6），与上面的偏导数最多差 1.3×10−10；取线性场 u=0.3+0.7x−0.4y（壁面值也取自它），修正后的扩散通量与 −νn·∇u 最多差 5.6×10−16（本页代码在构建时计算）。
Check: on the 40 interior and wall faces of this page, central differences of faceFlux() (step 10−6) differ from the partial derivatives above by at most 1.3×10−10; for the linear field u=0.3+0.7x−0.4y (the wall values taken from it as well), the corrected diffusive flux differs from −νn·∇u by at most 5.6×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 f17 (L=c9, R=c10). 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) 与扩散的两点部分形式相同，相当于数值扩散系数 ε|n|/p；在这张粗网格的内部面上，它介于 0.392 与 0.648 之间（壁面面上没有耗散项），平均每加密一次约减半。在收敛解上，ν=0.5 时各内部面耗散通量的绝对值之和是扩散通量的 0.87 倍；若取 ν=0.1，这个比值是 4.6，物理扩散会被数值耗散淹没。
Why ν=0.5. The dissipation −ε|n|(uR−uL) has the form of the two-point part of diffusion, a numerical diffusivity ε|n|/p; on the interior faces of this coarse mesh it lies between 0.392 and 0.648 (the wall faces carry no dissipation), and on average it roughly halves with every refinement. At the converged solution with ν=0.5, the absolute dissipative fluxes over the interior faces sum to 0.87 times the diffusive ones; with ν=0.1 the ratio would be 4.6, and the physical diffusion would be swamped by the numerical one.
到这里，内部面上的通量都有了。边界面还差一步：壁面面外侧没有单元，右状态只能是给定的壁面值 uw；远场面外侧也没有单元，外侧状态要由来流值或外推补上。第 4 节把这两种面补齐。
By now every interior face has its flux. The boundary faces need one more step: a wall face has no cell beyond it, so its right state can only be the given wall value uw; a far-field face has no cell beyond it either, and its outer state comes from the free stream or from extrapolation. Section 4 completes both kinds of face.
物面上 u=uw=0。格心格式在物面上没有未知量，壁面值只能从壁面面进入壁面单元的守恒方程（图 4.1）。壁面面上的状态已知，右状态就是 uw，通量是
On the body u=uw=0. A cell-centred scheme has no unknown on the body, so the wall value can only enter the wall cell's conservation equation through the wall face (Figure 4.1). The state on a wall face is known, the right state is uw itself, and the flux is
对流部分直接用已知的壁面状态：12uw2vn，没有平均，也没有耗散；这里 uw=0，所以它为零，穿过物面的只有扩散。扩散部分还是同一个修正公式，只是 ξ 从壁面单元的顶点平均点 xL 指向壁面面的中点，而平均梯度就是壁面单元自己的梯度 GL。此外，壁面值还作为一个模板点进入 GL。壁面就从这两处进入方程，都在普通的循环里：没有哪一行是特殊的，24 个方程全是守恒方程。
The convective part uses the known wall state directly, 12uw2vn, with no average and no dissipation; here uw=0 makes it zero, and only diffusion crosses the body. The diffusive part is the same corrected formula, except that ξ runs from the wall cell's vertex average xL to the wall face's midpoint and the averaged gradient is the wall cell's own GL. Besides, the wall value enters GL as a stencil point. Those are the wall's two doors into the equations, both inside the ordinary loops: no row is special, and all 24 equations are conservation equations.
与格点页面对照：那里壁面节点上 ui=uw 严格成立，但壁面通量要从被替换的行里找回来；这里贴壁单元的值不等于 uw，壁面值是经由通量“拉”过去的，只在网格加密的极限下成立。反过来，壁面通量本身就是残差里的一个普通面通量，直接就是目标函数（第 7 节）。
Compared with the node-centred page: there ui=uw holds exactly at the wall nodes, but the wall flux has to be recovered from the replaced rows; here the values of the wall cells are not uw: the wall value is pulled in through the flux, and u=uw holds only in the limit of refinement. In return, the wall flux is an ordinary face flux of the residual, and it is directly the objective (Section 7).
远场面上只算对流通量，扩散通量取零。对流部分与对流页面相同：远场面按 vn 的符号分成入流和出流（在初始几何上冻结），入流取右状态 uR=u∞=1，出流取 uR=uL。8 个远场面中，4 个入流，4 个出流。图 4.2 把每个远场面上的通量拆开：中心对流部分 14vn(uL2+uR2)，耗散部分 −ε|n|(uR−uL)，扩散部分 0。
The far-field 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 faces are split into inflow and outflow by the sign of vn (frozen at the initial geometry), the right state being uR=u∞=1 on inflow and uR=uL on outflow. Of the 8 far-field faces, 4 are inflow and 4 outflow. Figure 4.2 takes the flux on every far-field 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 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).
本页的解正是这样。在本页网格上，4 个入流面所在的单元上，u 介于 0.734 与 0.879 之间，而不是 1。一致加密（第 15 节的网格族）到第 5 层，入流面所在的单元上的 u 仍介于 0.67 与 0.97 之间；由最小二乘梯度量出的外法向导数 ∂nu 介于 0.06 与 0.56 之间，并不趋于零，而且与上式由 u 算出的值越来越近，两者之差在各个入流单元上的均方根逐层是 0.144、0.064、0.024、0.009、0.005、0.003（逐点的最大差在最细两层是 0.014、0.015，没有再减小；它出现在入流转为出流的那个角点旁）。出流一侧的 ∂nu 则趋于零：它的最大绝对值从第 1 层起逐层是 0.103、0.080、0.053、0.035、0.025（图 4.3）。附录 A.1 第 3 小节正是从这个总通量条件出发，推导远场的伴随边界条件。
The page's solution does exactly this. On the page's mesh, in the cells of the 4 inflow faces, u lies between 0.734 and 0.879, not at 1. Refined uniformly (the family of Section 15) to level 5, u in the cells of the inflow faces still lies between 0.67 and 0.97; the outward normal derivative ∂nu measured by the least-squares gradient lies between 0.06 and 0.56, 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 cells being 0.144, 0.064, 0.024, 0.009, 0.005 and 0.003 level by level (the largest pointwise difference, 0.014 and 0.015 on the two finest levels, no longer falls; it sits next to the corner where inflow turns into outflow). On the outflow side ∂nu does go to zero: its largest magnitude is 0.103, 0.080, 0.053, 0.035 and 0.025 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.
为什么对流页面没有这种壁面。纯对流方程是一阶的，边界上只能在特征进入区域的地方给数据，而特征速度沿闭合的物面必然有正有负（见对流页面第 4 节）。扩散项是二阶的，它让整条物面都可以、也必须给一个条件。
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, and along a closed body the characteristic speed takes both signs (see Section 4 of the advection page). The diffusion term is second order: it lets, and requires, the whole body take one condition.
本页每个单元存储 1 个未知量。以下用 表示把所有位置按代码顺序排成的列向量（标量问题每个位置只有一个分量）。它是一个有限维向量，不是连续函数。物理方程中的空间导数只用于说明模型；求离散伴随时，求导对象是本页实际装配的有限维函数。
Each cell 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 cell, the wall value of each wall face 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 diffusive flux needs the gradients of the cells on either side of a face, and a cell gradient needs the neighbours' values, so the gradients of all cells must be finished before the first flux is computed. The residual therefore becomes two face loops:
没有第三步：弱施加不替换任何一行。图 5.1 把三个函数画成数据流：前向走同样的顺序，伴随把顺序倒过来。
There is no third step: weak imposition replaces no row. Figure 5.1 draws the three functions as data flows: the forward code runs in the same order, the adjoint reverses it.
G_b 里，再把梯度那一遍反着走；uw 的那一份丢弃（对壁值求导时就收集它，见第 11 节）。Figure 5.1 · The data flow of the three functions. Top: residual, which computes the gradients G from u (and the wall value uw), then r from u and G. Middle: residual_d, the same order acting on the tangent v; uw is a constant, whose tangent is zero. Bottom: residual_b, in reverse: the flux pass runs backwards first, parking the gradients' share in the array G_b, then the gradient pass runs backwards; the share of uw is dropped (a derivative with respect to the wall value collects it; see Section 11).图 5.2 逐句执行这两遍循环。单元里的紫色箭头是梯度数组 G：第一遍里一个面一个面地长出来，壁面面上的差值来自面中点的 uw；第二遍里只读不写。下方三行数组依次是状态 u*、梯度 G 和输出 r。
Figure 5.2 executes both loops statement by statement. The violet arrows in the cells are the gradient array G: in pass 1 it grows face by face, the wall faces taking their difference from the uw at the face midpoint; in pass 2 it is only read. The three arrays below are the state u*, the gradients G and the output r.
播放到第二遍的某个内部面时，留意它收集了两侧的状态和两侧的梯度，而每个梯度又来自整个模板，所以一个面的通量依赖两个模板里的所有单元。播放到壁面面时，右状态是面中点上的 uw，平均梯度只有 GL。循环结束时 r 在收敛解上处处为零，包括 8 个壁面单元：它们的方程和其他单元一样，是一个完整的守恒方程。
When the player reaches an interior face in pass 2, note that it gathers both states and both gradients, and each gradient came from a whole stencil, so one face's flux depends on every cell of the two stencils. On a wall face the right state is the uw at the face midpoint and the averaged gradient is GL alone. When the loop ends, r is zero everywhere at the converged solution, the 8 wall cells included: their equations, like every other cell's, are complete conservation equations.
这一节先推出残差的精确导数：切线方程与伴随方程的矩阵正是它（第 8、9 节不形成矩阵，直接作用它）；再用 Jacobi 迭代求解 R(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 R(u)=0 with the Jacobi iteration, which needs only its diagonal.
对流部分只依赖面两侧的单元值：∂ha/∂uL=12vnuL+ε|n|=cL，∂ha/∂uR=12vnuR−ε|n|=cR；壁面面上 hw 的对流部分与 uL 无关，cL=0。扩散部分有两条路径：两点差直接依赖 uL、uR，平均梯度经过最小二乘依赖两个模板中的每一个单元（下式中 [m=R] 在 m=R 时为 1，否则为 0）：
The convective part depends only on the cell values either side of the face: ∂ha/∂uL=12vnuL+ε|n|=cL and ∂ha/∂uR=12vnuR−ε|n|=cR; on a wall face the convective part of hw does not depend on uL, so cL=0. The diffusive part has two paths: the two-point difference depends directly on uL and uR, the averaged gradient goes through least squares and depends on every cell of the two stencils (below, [m=R] is 1 if m=R and 0 otherwise):
（壁面面上没有 R，梯度只有 GL，系数 12 换成 1；壁面中点那个模板点是常数，对它没有导数。）所以精确 Jacobian 的每一行伸到面邻居的面邻居：壁面单元 10 个非零元，三角形单元 7 或 8 个，全矩阵 200 个。
(A wall face has no R, its gradient is GL alone and the 12 becomes 1; the wall-midpoint stencil point is a constant, with no derivative.) So each row of the exact Jacobian reaches the face neighbours of the face neighbours: 10 non-zeros for a wall cell, 7 or 8 for a triangle, 200 in the whole matrix.
只保留相邻单元之间的耦合，就得到精确 Jacobian 的一阶部分 A1：对流导数原样保留，扩散导数只保留两点部分，梯度路径整个丢掉。很多隐式求解器就拿它作迭代矩阵。本页的求解不用它：它在这里只用来看清梯度路径添了什么（图 6.1），以及第 8 节里，拿它的转置求伴随会错多少。
Keep only the couplings between adjacent cells and what is left is the first-order part A1 of the exact Jacobian: the convective derivatives as they are, and of diffusion only the two-point part, the gradient path dropped. 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 的非零结构：每行 3 或 4 个，全矩阵 88 个。
A1 has exactly the non-zero pattern of the advection page's Jacobian: 3 or 4 per row, 88 in all.
第一步，一个面的通量读到哪些未知量。内部面的通量读两个状态 uL、uR 和两个梯度 GL、GR；每个梯度又读它自己的单元和它模板里的单元。所以 hf 依赖 L、R 以及两个模板里的所有单元（图 6.2 左，面 f17：5 个）。壁面面只有左侧：cL 和它的模板（右，4 个，外加数据 uw）。远场面没有扩散通量，只读 uL。
Step 1, what one face flux reads. The flux of an interior face reads two states, uL and uR, and two gradients, GL and GR; each gradient in turn reads its own cell and the cells of its stencil. So hf depends on L, R and every cell of the two stencils (Figure 6.2, left: the face f17, 5 cells). A wall face has a left side only: cL and its stencil (right, 4 cells, plus the data uw). A far-field face carries no diffusive flux and reads uL only.
第二步，一行是它所有面的并集。第 c 行收集 ±∂h/∂um，求和取遍单元 c 的所有面，所以它的非零列就是这些面读到的未知量的并集：c 自己、它的面邻居，以及邻居的模板成员，也就是邻居的邻居。
Step 2, a row is the union over its faces. Row c collects ±∂h/∂um over all faces of cell c, so its non-zero columns are the union of what those faces read: c itself, its face neighbours, and the members of the neighbours' stencils, the neighbours of neighbours.
第三步，按单元的种类数。三种单元各 8 个，同一种单元的行非零元个数相同：壁面单元 10 个，内部三角形 8 个，外圈三角形 7 个（它的远场面不进模板，邻居少一个）。一阶部分 A1 只保留自己和面邻居：4、4、3 个。
Step 3, count by kind of cell. There are 8 cells of each kind, and the rows of one kind have the same number of non-zeros: 10 for a wall cell, 8 for an inner triangle, 7 for an outer triangle (its far-field face is in no stencil, so it has one neighbour fewer). The first-order part A1 keeps only the cell and its face neighbours: 4, 4 and 3.
核对：按第一、二步从每个面读到的集合拼出的 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），第一遍只点亮碰到它的面，把差值分给两侧的梯度，第二遍再点亮所有碰到这个状态或被改变梯度的面，写进它们两侧的行，最后得到 Av=Aêj，即 A 的第 j 列。反向：给第 i 个残差一个伴随种子（w=êi），第二遍倒过来只读第 i 行的面，第一遍倒过来再经最小二乘权重读到模板里的每个单元，最后得到 ATêi，即 A 的第 i 行。默认选 c10，它的行正是 图 6.1 左边那一行。本页上同一个单元的行与列非零位置相同（c10 行 8 个、列 8 个），数值不同：A 不对称，所以伴随要一个与前向不同的函数。
Figure 6.3 sets this structure in motion. Forward: seed one cell with a tangent (v=êj); pass 1 lights up only the faces that touch it and hands the difference to the gradients on both sides, pass 2 then lights up every face that touches the seeded state or a changed gradient and writes into the rows on both of its sides, and the result is Av=Aêj, column j of A. Reverse: give residual i an adjoint seed (w=êi); pass 2 reversed reads only the faces of row i, pass 1 reversed then reaches, through the least-squares weights, every cell of the stencils, and the result is ATêi, row i of A. The default is c10, whose row is exactly the one on the left of Figure 6.1. On this page the row and the column of one cell have the same non-zero positions (8 in row 10, 8 in column 10) but different values: A is not symmetric, which is why the adjoint needs a function of its own.
与格点页面一样，本页的 primal、切线、伴随三个方程（后两个在第 10、11 节）用的是同一个 Jacobi 迭代：每一轮算一次残差，再逐单元除以精确 Jacobian 的对角元：
As on the node-centred 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, cell by cell, by the diagonal of the exact Jacobian:
diag(A) 不需要把 A 造出来：diagonal(u) 就是残差的那个面循环，每个面只留下落在对角线上的两项，即对流导数、两点部分，以及经由最小二乘梯度的那一项（壁面面上系数 12 换成 1，只经由 GL）。所以一轮只有一次 residual 和一次 diagonal：不形成矩阵，也不解线性方程组。没有替换行：壁面是弱施加的，每个单元都照常更新。初值取 u≡1，判据是残差的均方根 ‖R‖<10−13（图 6.4）。
diag(A) does not need A to be formed: diagonal(u) is the residual's own face loop, keeping of each 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 (on a wall face with 1 in place of 12, and through GL only). So a sweep is one residual and one diagonal: no matrix is formed and no linear system is solved. There are no replaced rows: the wall is weak, and every cell is updated alike. The initial guess is u≡1, and the criterion is the rms of the residual, ‖R‖<10−13 (Figure 6.4).
为什么用 Jacobi，而不用 LU。本页只有 24 个未知量，理论上最快的办法是直接分解：把精确的 A 整个造出来，用 Newton 法，每一步做一次稠密 LU 分解，k=5 时残差就已低于 10−13。可那样就说明不了 matrix-free。Jacobi 是最简单的迭代法：对矩阵它只要对角元，其余的一切都经由残差作用进来（切线方程里是 residual_d，伴随方程里是 residual_b），这正是 matrix-free 的意思。代价是线性收敛：在解附近，误差每一轮乘以迭代矩阵 I−diag(A)−1A，速率由它的谱半径决定，本算例是 0.915，判据在 k=315 时满足。切线与伴随的迭代矩阵与它同谱，所以进入渐近阶段以后，三条收敛曲线以同一个速率下降（第 11 节）。
Why Jacobi, and not LU. With only 24 unknowns the fastest method in principle is a direct factorisation: form the whole exact A, run Newton, and factorise with one dense LU per step; the residual is below 10−13 by k=5. 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(A)−1A, at a rate set by its spectral radius, 0.915 here, and the criterion is met at k=315. 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 的对角元，迭代照样收敛到同一个 u*（残差是精确的），只是速率变了。在这里差别很大：梯度路径能改变一个对角元的 15.3%，谱半径从 0.915 升到 0.990，要多跑九倍左右的轮次。格点页面上它最多只改变 3.9%，谱半径几乎不动：精确对角元时 0.787，取一阶部分的对角元时 0.786。
The diagonal must carry the gradient path. Take the diagonal of the first-order part A1 instead and the iteration still converges to the same u* (the residual is exact), only at another rate. Here the difference is large: the gradient path changes a diagonal entry by up to 15.3%, and the spectral radius rises from 0.915 to 0.990, about nine times as many sweeps. On the node-centred page it changes one by at most 3.9%, and the spectral radius hardly moves: 0.787 with the exact diagonal, 0.786 with the first-order one.
固定轮次，但会检查。代码固定跑 600 轮，而不是一满足判据就退出：余下的轮次把 u* 一直带到舍入误差平台（从 k=386 起），页面上每一个实测常数都是在这些轮次下量出来的；得到的 u* 与 Newton 法的解最多差 3.3×10−16。它记录判据在第几轮满足，若 600 轮内从未满足就直接报错，而不会把一个未收敛的状态当成 u* 返回。
A fixed sweep count, but a checked one. The code runs a fixed 600 sweeps rather than exiting at the criterion: the remaining sweeps take u* down to its round-off floor (from k=386 on), and every measured constant on this page was produced under those sweeps; the u* it returns is within 3.3×10−16 of Newton's. It records the sweep at which the criterion was met, and raises an error if 600 sweeps never meet it, rather than returning an unconverged state as u*.
| 初值处的 ‖R‖ | ‖R‖ at the initial guess | 0.489 |
|---|---|---|
| 谱半径 ρ(I − diag(A)−1A) | spectral radius ρ(I − diag(A)−1A) | 0.915 |
| 满足 ‖R‖ < 10⁻¹³ 的轮次 | sweep meeting ‖R‖ < 10⁻¹³ | k = 315（Newton 加稠密 LU：5）k = 315 (Newton with a dense LU: 5) |
| 舍入误差平台的起点 | where the round-off floor starts | k = 386 |
| 600 轮之后的 ‖R‖ | ‖R‖ after the 600 sweeps | 2.7e−16 |
| u* 的范围 | the range of u* | 0.179 – 0.879 |
第一步，把一轮写成一个映射：uk+1=uk−diag(A(uk))−1R(uk)。diag(A) 随 u 变，因为精确 Jacobian 随 u 变（对流导数 cL、cR 里有 u）。
Step 1, write one sweep as a map: uk+1=uk−diag(A(uk))−1R(uk). diag(A) changes with u, because the exact Jacobian does (the convective derivatives cL, cR contain u).
第二步，在收敛解附近线性化。记第 k 轮的误差为 ek=uk−u*，并用 R(u*)=0：
Step 2, linearise near the converged solution. Write ek=uk−u* for the error after k sweeps and use R(u*)=0:
diag(A) 的变化乘的是 R(u*)=0，所以在一次项里消失了：收敛解附近，一轮就是把误差乘以迭代矩阵 I−diag(A)−1A。
The change of diag(A) multiplies R(u*)=0 and drops out of the linear term: near the converged solution one sweep multiplies the error by the iteration matrix I−diag(A)−1A.
第三步，重复 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.
第四步，一个充分条件，以及它为什么在这里不成立。迭代矩阵的对角元为零，非对角元是 −Acj/Acc。Gershgorin 圆盘定理说，每个特征值都落在某个以 0 为心、半径为 σc 的圆盘里：
Step 4, a sufficient condition, and why it fails here. The iteration matrix has a zero diagonal and off-diagonal entries −Acj/Acc. The Gershgorin disc theorem puts every eigenvalue inside one of the discs around 0 of radius σc:
若每一行都严格对角占优（σc<1），收敛就有保证。本页却有 12 行 σc>1，最大 1.154：这个界说明不了什么。收敛要看特征值本身（图 6.5）：它们全在单位圆内，ρ=0.9149，第二大的模是 0.8561。
If every row were strictly diagonally dominant (σc<1), convergence would be guaranteed. Here 12 rows have σc>1, the largest 1.154: the bound says nothing. Convergence rests on the eigenvalues themselves (Figure 6.5): all of them lie inside the unit circle, ρ=0.9149, and the second largest modulus is 0.8561.
第五步，对角元里的梯度路径。改用一阶部分的对角元 diag(A1)，迭代矩阵变成 I−diag(A1)−1A。梯度那一项在外圈三角形上最大（图 6.6，最多占对角元的 15.3%），而且多为正：diag(A1) 比 diag(A) 小，除得少，步子就迈大了。于是出现一个接近 −1 的特征值 −0.9905：误差每一轮几乎原样翻号，来回振荡，谱半径升到 0.9905。
Step 5, the gradient path in the diagonal. With the first-order part's diagonal diag(A1) the iteration matrix becomes I−diag(A1)−1A. The gradient term is largest on the outer triangles (Figure 6.6, up to 15.3% of a diagonal entry) and mostly positive: diag(A1) is smaller than diag(A), divides by less, and so takes steps that are too long. The result is an eigenvalue close to −1, −0.9905: the error flips its sign almost unchanged every sweep, oscillating, and the spectral radius rises to 0.9905.
核对：在 u* 附近对一轮 Jacobi 映射取中心差分（步长 10−6），与 I−diag(A)−1A 最多差 5.6×10−11；由特征值算出的 ρ=0.914894，迭代在第 150 到 250 轮实测的收缩因子是 0.914894（本页代码在构建时计算）。
Check: central differences of one Jacobi sweep near u* (step 10−6) differ from I−diag(A)−1A by at most 5.6×10−11; the eigenvalues give ρ=0.914894, and the iteration shrinks by 0.914894 per sweep measured over k=150 to 250 (computed by this page's code at build time).
弱施加让目标函数变得直接：物面吸收的通量就是 8 个壁面面（集合 Fw）上的通量，它们本来就在残差里：
Weak imposition makes the objective direct: the flux the body absorbs is the flux through the 8 wall faces (the set Fw), which are part of the residual already:
它自动与离散格式一致：把 24 个单元的守恒方程加起来，内部面的通量两两抵消，剩下的就是壁面面与远场面通量之和，所以 J 恰好等于从远场流入的净通量。实测 J=2.343758，远场面通量之和是 −2.343758，两者相加为 0。格点页面要从被替换的行里把这个量找回来，这里不需要。
It is consistent with the discrete scheme automatically: add up the conservation equations of all 24 cells, the interior faces cancel in pairs, and what is left is the sum over the wall and the far-field faces, so J is exactly the net flux entering through the far field. Measured: J=2.343758, the far-field faces sum to −2.343758, and the two add up to 0. The node-centred page has to recover this quantity from its replaced rows; here there is nothing to recover.
每个壁面面的通量里，两点部分 νp(uL−uw) 占绝大部分，梯度修正 −νGL·k 很小：只取两点部分，总和是 2.305925，修正贡献了 1.6%（图 7.1）。
In each wall face's flux the two-point part νp(uL−uw) carries almost everything and the gradient correction −νGL·k is small: the two-point parts alone sum to 2.305925, the correction adding 1.6% (Figure 7.1).
两种格式在这张粗网格上给出的壁面通量差别不小：格心 2.344，格点 2.749。一致加密（第 15 节的网格族）之后，格心依次为 2.344、2.331、2.249、2.184、2.145、2.124（24、96、384、1536、6144、24576 个单元），格点在同一族网格上依次为 2.749、2.454、2.291、2.201、2.153、2.127：两者向同一个值靠拢。
On this coarse mesh the two schemes give quite different wall fluxes: 2.344 cell-centred, 2.749 node-centred. Refine uniformly (the family of Section 15) and the cell-centred value runs 2.344, 2.331, 2.249, 2.184, 2.145 and 2.124 (24, 96, 384, 1536, 6144 and 24576 cells), while the node-centred one on the same family runs 2.749, 2.454, 2.291, 2.201, 2.153 and 2.127: the two close in on the same value.
J 是 8 个面通量之和，所以 g=(∂J/∂u)T 就是一次伴随循环：壁面面的种子取 1、其余取 0，走一遍 residual_b 的面循环（第 9 节）。壁面面通量依赖壁面单元的值与它的梯度，梯度又依赖整个模板，所以 g 在 16 个单元上非零：8 个壁面单元，以及它们模板里的单元。
J is a sum of 8 face fluxes, so g=(∂J/∂u)T is one adjoint loop: seed 1 on the wall faces and 0 elsewhere, through the face loop of residual_b (Section 9). A wall-face flux depends on the wall cell's value and on its gradient, and the gradient on the whole stencil, so g is non-zero on 16 cells: the 8 wall cells and the cells in their stencils.
第一步，目标函数是 8 个壁面面通量之和，每一项读壁面单元的值和它的梯度：
Step 1, the objective is the sum of the 8 wall-face fluxes, each reading the wall cell's value and its gradient:
第二步，逐面求导，再经最小二乘梯度链式求导（第 3 节两段推导的第三、五步）：
Step 2, differentiate face by face, and chain through the least-squares gradient (steps 3 and 5 of the two derivations of Section 3):
第三步，读出各分量。m 是壁面面 f 自己的单元 Lf 时，∂GL/∂uL=−ΣswLs；m 是 Lf 模板里的邻居时，∂GL/∂um=wLm；其余为零：
Step 3, read off the components. For m the cell Lf of the wall face f itself, ∂GL/∂uL=−ΣswLs; for m a neighbour in the stencil of Lf, ∂GL/∂um=wLm; otherwise zero:
所以 g 只在壁面单元和它们模板里的单元上非零，共 16 个（图 7.2）。壁面单元上的大数主要来自两点部分 νpf；模板单元上的小数来自经由 GL 的修正。
So g is non-zero only on the wall cells and the cells in their stencils, 16 in all (Figure 7.2). The large values in the wall cells come mostly from the two-point part νpf; the small ones in the stencil cells from the correction through GL.
第四步，代码怎么算。上式正是 residual_b 的两遍反向扫描，面通量的伴随取壁面面上 h̄=1、其余面上 0：第二遍倒过来给 ūL 加上 νp、给 ḠL 加上 −νk，第一遍倒过来把 Ḡ 经权重送回模板。
Step 4, how the code computes it. The formula above is exactly the two reverse passes of residual_b with the face-flux adjoint h̄=1 on the wall faces and 0 elsewhere: pass 2 reversed adds νp to ūL and −νk to ḠL, and pass 1 reversed sends Ḡ back to the stencil through the weights.
第五步，对坐标的偏导（状态不动）。同一个 J=Σf∈Fwhw 也直接依赖坐标（n、ξ、p、k 与最小二乘权重）：
Step 5, the partial derivative with respect to the coordinates, the state held fixed. The same J=Σf∈Fwhw depends on the coordinates directly too (through n, ξ, p, k and the least-squares weights):
它是一次对坐标的反向扫描，壁面面上的种子为 1：第 12 节“壁面面再加 1”说的就是这一项。
It is one reverse sweep over the coordinates with the seed 1 on the wall faces: the “1 more on every wall face” of Section 12 is exactly this term.
核对：按上面的公式逐项算出的 g 与代码 dJdu()（即 residual_b）逐位相同，与对 J 的中心差分（步长 10−6）最多差 2.2×10−10；∂J/∂X 的反向扫描与状态不动时移动坐标的中心差分最多差 3.4×10−10，它的分量最大为 1.295（本页代码在构建时计算）。
Check: the g computed term by term from the formulas above agrees with the code's dJdu() (that is, residual_b) to the last bit and differs from central differences of J (step 10−6) by at most 2.2×10−10; the reverse sweep for ∂J/∂X differs from central differences that move the coordinates with the state held fixed by at most 3.4×10−10, whose largest component is 1.295 (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 of this page is such a weighted sum: one term per wall face, and ; a wall-face flux depends on the wall cell's value, its least-squares gradient, the wall value and the geometry, so depends on both the state and the design variables.
求 时固定所有设计变量；求 时固定全部离散状态槽位。固定解析目标函数并不一定意味着移动网格上的目标采样值不变：如果采样点随设计移动，还要保留该位置依赖；固定网格及固定采样值时，相应项才为零。目标中的体积权重已经进入 ，不应在伴随右端再乘一次体积。
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 节末尾逐类给出。对每一组 X，收敛解都满足 R(u(X),X)=0，所以它对 X 的导数也是零；再用链式法则写出目标函数的全导数（g=(∂J/∂u)T，第 7 节）：
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. For every X the converged solution satisfies R(u(X),X)=0, so its derivative with respect to X is zero as well; and the chain rule gives the total derivative of the objective (g=(∂J/∂u)T, Section 7):
du/dX 是 24×48 的矩阵，直接算它要解 48 次线性方程组，每个坐标一次（第 10 节的切线方程）。把乘积换一个结合顺序，先算行向量 gTA−1：它与设计变量无关，只要解一次方程。记 ψT=−gTA−1，也就是 ATψ=−g，就得到
du/dX is a 24×48 matrix, and computing it directly takes 48 linear solves, one per coordinate (the tangent equation of Section 10). Associate the product the other way and compute the row vector gTA−1 first: it does not depend on the design variables and costs one solve. Write ψT=−gTA−1, that is ATψ=−g, and
这就是第 11、12 节做的事：一次伴随求解，加一次对坐标的反向扫描。同一个推导用有限维 Lagrangian 写出来，在第 11 节末尾。
That is what Sections 11 and 12 do: one adjoint solve and one reverse sweep over the coordinates. The same derivation, written with a finite-dimensional Lagrangian, is at the end of Section 11.
伴随方程里的矩阵是精确的 A，不是它的一阶部分 A1（第 6 节；很多隐式求解器留用的就是它）。拿一阶部分的转置 A1T 解出来的 ψ 满足的是另一个方程：本算例里它与精确的 ψ 最多差 0.040，给出的形状梯度最多差 1.6e−2（最大分量的 3.8%），而精确伴随与有限差分的差在 10−9 量级（第 13 节）。精确矩阵每行伸到邻居的邻居，不形成它，只需要能算 Av 与 ATw 的两个函数：residual_d 与 residual_b。
The matrix of the adjoint equation is the exact A, not its first-order part A1 (Section 6), which many implicit solvers keep. A ψ solved with the transpose of the first-order part, A1T, satisfies a different equation: here it misses the exact ψ by up to 0.040, and its shape gradient misses by up to 1.6e−2 (3.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; it is not formed, and all that is needed are the two functions that compute Av and ATw: residual_d and residual_b.
记号。本页的前向与伴随都是对离散代码求导，用自动微分的记号：程序变量 x 的切线（前向模式的导数）记作 ẋ，伴随（反向模式的导数）记作 x̄；伴随先清零（x̄←0），再按语句的逆序累加（x̄+=…）。代码与演示里用 Tapenade 的名字 xd≡ẋ、xb≡x̄：residual_d 的输入 v 就是状态的切线 u̇，输出 r_d 就是 ṙ=Au̇。δ 只表示连续函数的变分，只有附录 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 ṙ=Au̇. The symbol δ is kept for the variation of a continuous function: only the continuous adjoint of Appendix A uses the variation symbol.
第一遍的梯度对 u 线性，它的切线就是同一个循环作用在 u̇ 上，只是壁面中点上的 uw 是常数，切线为 0。第二遍按两侧状态与平均梯度 Gf 求导；壁面面上 cL=0，u̇R=0。没有替换行，也就没有“循环之后覆盖”这一步：
The gradients of pass 1 are linear in u, so their tangent is the same loop acting on u̇, except that the uw at a wall midpoint is a constant, whose tangent is zero. Pass 2 is differentiated with respect to the two states and the averaged gradient Gf; on a wall face cL=0 and u̇R=0. There are no replaced rows, and so no "overwrite after the loop":
每条从 u 到 r 的路径都被求了导，这就是精确的 ṙ=Au̇，即 Av。去掉 Ġf·k，剩下的就是一阶部分 A1v。图 8.1 逐句执行 residual_d，并与显式的 A（只为校对而形成）对照。
Every path from u to r is differentiated: this is the exact ṙ=Au̇, that is Av. Drop Ġf·k and what is left is the first-order part A1v. Figure 8.1 executes residual_d statement by statement, against the explicit A (formed for the check only).
由链式法则，residual_d(u,v) 就是 ddtR(u+tv) 在 t=0 处的值。播放到最后，矩阵视图下方给出它与显式乘积的最大差，在舍入量级。
By the chain rule, residual_d(u,v) is ddtR(u+tv) at t=0. Play to the end and the matrix view reports its largest difference from the explicit product, at round-off level.
上面的切线公式已经够算 Au̇。第 10 节的右端项和第 12 节的设计导数还要更多：每一种边界面，以及设计输入（来流值、源项、壁值、坐标）进入残差的每一处。下面把 residual 的每一句（第 5 节的伪代码）单独求导。只用一条规则，就是链式法则：一句的每个输入 x 带着自己的切线 ẋ，输出的切线是各输入的切线乘这一句的偏导数再相加。
The tangent formulas above are enough for Au̇. 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 行）：差值 Δ=uR−uL（壁面面上 uR=uw），再累加 GL+=wLRΔ 与 GR+=−wRLΔ。输入是两个状态和两个权重，所以
Step 1, pass 1 (lines 4 to 6): the difference Δ=uR−uL (uR=uw on a wall face), and the two accumulations GL+=wLRΔ and GR+=−wRLΔ. The inputs are two states and two weights, so
权重只依赖几何，ẇ 只在坐标带切线时不为零（第 12 节推出它）。壁面面上 u̇R=u̇w，对状态求切线时它是 0。只对状态求切线时只剩各式的第一项：第一遍的切线就是同一个循环作用在 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). On a wall face u̇R=u̇w, zero in the state's tangent. 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.
第二步，第二遍的收集（第 9、10 行）：取两侧状态与平均梯度，切线照样取。右状态随面的种类而定：内部面上 u̇R 是右侧单元的切线，壁面面上 u̇R=u̇w，远场入流面上 u̇R=u̇∞，出流面上 uR=uL，所以 u̇R=u̇L。平均梯度 Ġf=12(ĠL+ĠR)，壁面面上 Ġf=ĠL。对状态求切线时 uw 与 u∞ 是常数，切线为 0。
Step 2, the gathers of pass 2 (lines 9 and 10): the two states and the averaged gradient, and their tangents likewise. The right state depends on the kind of face: on an interior face u̇R is the tangent of the cell on the right, on a wall face u̇R=u̇w, on a far-field inflow face u̇R=u̇∞, and on an outflow face uR=uL, so there u̇R=u̇L. The averaged gradient gives Ġf=12(ĠL+ĠR), and Ġf=ĠL on a wall face. In the state's tangent uw and u∞ are constants, with zero tangents.
第三步，对流通量（第 11 行的前一半）。vn=β·n 与 |n| 都随 n 变：v̇n=β·ṅ，ddt|n|=n·ṅ/|n|。于是
Step 3, the convective flux (the first half of line 11). Both vn=β·n and |n| move with n: v̇n=β·ṅ and ddt|n|=n·ṅ/|n|. Hence
远场面只是这条式子的特例；壁面面不是：那里状态已知，ha=12uw2vn，没有耗散项，ḣa=uwvnu̇w+12uw2β·ṅ，与 uL 无关（cL=0，cR=uwvn）。出流面上 uR=uL（第 4 节），两个槽位合成 (cL+cR)u̇L=vnuLu̇L：代码在出流面上的 cL 就是这个合并后的系数。入流面上的 cRu̇∞ 是 u∞ 进入残差的唯一一处，它就是第 10 节的 f。
The far-field faces are special cases of this formula; the wall face is not: there the state is known, ha=12uw2vn with no dissipation, and ḣa=uwvnu̇w+12uw2β·ṅ, independent of uL (cL=0, cR=uwvn). On an outflow face uR=uL (Section 4), and the two slots merge into (cL+cR)u̇L=vnuLu̇L: the code's cL on an outflow face is this merged coefficient. On an inflow face the term cRu̇∞ is the one place where u∞ enters the residual: it is the f of Section 10.
第四步，扩散通量（第 11 行的后一半）。hv 读四样东西：两侧状态、平均梯度、p 和 k；而 p=ξ·n/|ξ|2、k=n−pξ 由 n 与 ξ 算出（第 3 节）：
Step 4, the diffusive flux (the second half of line 11). 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∂ξ·ξ̇。壁面面上 uR=uw、Gf=GL，状态部分多出 −νpu̇w。远场面上 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 wall face uR=uw and Gf=GL, and the state part gains −νpu̇w. On a far-field face hv=0, and so is its tangent.
第五步，分发（第 12、13 行）与源项。累加是线性的，切线照样累加；源项 sc（第 11 节的设计变量）进入 Rc←Rc−sc：
Step 5, the scatters (lines 12 and 13) and the source. Accumulations are linear, and their tangents accumulate the same way; a source sc (a design variable of Section 11) enters as Rc←Rc−sc:
第六步，哪个种子喂哪一句。把第一至五步按输入读一遍，就知道每一类设计分量的右端项 f=∂R/∂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=∂R/∂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 满足 Av=−f，所以跑完以后只有 4 个入流单元上 ṙ=−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 face by face, pass 2 the fluxes' tangents face by face, and the bars on the right are the output ṙ accumulated in each cell. Since v satisfies Av=−f, at the end ṙ=−f is non-zero only on the 4 inflow cells 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.4×10−16，最小二乘权重的切线最多差 8.9×10−16；壁值种子下与中心差分（步长 10−4）最多差 5.8×10−13；源项种子与中心差分逐位相同；动画最后一帧与 −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.4×10−16 relative, and the least-squares weights' tangents by at most 8.9×10−16; with a wall-value seed they differ from central differences (step 10−4) by at most 5.8×10−13; 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(A)−1(ATψ+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(A)−1(ATψ+g) can only stop where
答案由残差里的算子决定，与对角阵无关。用精确的 residual_b，但把 diag(A) 换成一阶部分的对角元 diag(A1)，迭代 6000 轮后与精确的 ψ 只差 4.2×10−15：第 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) in place of diag(A), 6000 sweeps end within 4.2×10−15 of the exact ψ: just as the primal of Section 6 reaches the same u* with another diagonal.
第二步，换掉算子，不动点就跟着换。若残差里用一阶部分，ψ̃ 解的是 A1Tψ̃=−g。减去 ATψ=−g，再在左边加减 A1Tψ：
Step 2, change the operator and the fixed point changes with it. With the first-order part in the residual, ψ̃ solves A1Tψ̃=−g. Subtract ATψ=−g and add and subtract A1Tψ on the left:
驱动误差的正是一阶部分丢掉的那一块，即梯度路径 A−A1；迭代多少轮、判据多严都去不掉它。
The error is driven by exactly what the first-order part drops, the gradient path A−A1; 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
核对：本页 ψ̃ 与 ψ 最多差 0.040，第二步的等式两边最多差 1.3×10−15；用 ψ̃ 算出的形状梯度最多差 1.6×10−2，是最大分量的 3.8%（图 8.3，本页代码在构建时计算）。
Check: on this page ψ̃ and ψ differ by at most 0.040, and the two sides of the equation of step 2 by at most 1.3×10−15; the shape gradient computed with ψ̃ misses by up to 1.6×10−2, 3.8% of the largest component (Figure 8.3; computed by this page's code at build time).
固定网格和当前可微分支，把一个面的所有输入记为 ；它包括重构状态以及该通量需要的梯度、几何和边界数据。令 、 分别是 对存储状态和设计变量的导数， 是局部通量 Jacobian。于是每个残差行的导数按同一个有向面关联矩阵 装配，这里的求和包括相应的边界通量映射。本页没有被替换的行，边界条件全部在边界面的通量里。前向模式按计算顺序传播切向量，状态种子是 ，设计种子是 ：
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. This page replaces no rows: every boundary condition sits in a boundary-face flux. 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 vertex averages on the cell 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, (cell 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.1786 | 0.2435 | |
| 0.2613 | 0.3573 | ||
| 0.2506 | 0.3501 | ||
| 0.3332 | 0.4916 | ||
| 0.3071 | 0.4068 | ||
| 设计种子design seed | 1 | 1 | |
| （2）重构(2) reconstruct | 0.2613 | 0.3573 | |
| 0.2506 | 0.3501 | ||
| −0.01891 −0.05842 | −0.03217 −0.06079 | ||
| 0.06801 0.146 | 0.08085 0.2311 | ||
| −1.366 0.0463 | 0 0 | ||
| 1.366 | 0 | ||
| −0.3195 0.2882 | 0 0 | ||
| （3）通量(3) flux | 0.002473 | −0.06985 | |
| （4）分发(4) scatter | 0.002473 | −0.06985 | |
| −0.002473 | 0.06985 |
不直接进入这个面：， 全部来自状态的切向量。 只作为远场入流面的右状态进入通量，所以 只在 4 个单元上非零；几何量也是 的一部分，但这个种子不移动任何节点，所以它们的切向量都是零。
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 4 cells 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 , which vanishes at the converged solution: the tangent equation solved further on.
按程序顺序执行四条语句，箭头上是各自的初等 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.2435， 0.3573， 0.3501， 0.4916， 0.4068。
Seeds: ( is the free-stream value ) and . On : 0.2435, 0.3573, 0.3501, 0.4916, 0.4068.
（1）取值：，只挑出 上的 5 个状态的切向量（上一步列出的数）。
(1) Gather: picks the tangents of the 5 states of (the numbers of the previous step).
（2）重构：。 0.3573， 0.3501， (−0.03217, −0.06079)， (0.08085, 0.2311)。 (0, 0)， 0， (0, 0)（这个种子不移动任何节点）。： 不进入这个面。
(2) Reconstruct: . 0.3573, 0.3501, (−0.03217, −0.06079), (0.08085, 0.2311). (0, 0), 0, (0, 0) (this seed moves no node). : does not enter this face.
（3）通量： −0.06985。
(3) Flux: −0.06985.
（4）分发： −0.06985，，即加上 0.06985。
(4) Scatter: −0.06985 and , that is 0.06985 added.
与格点页面相同，先把前向写成线性映射的乘积，再逐个转置。第一遍是 2N×N 的算子 Γ（N=24 是未知量个数）：Ġ=Γu̇，第 c 个块行是 Σswcs(ês−êc)T，其中壁面中点那一项只剩 −wcsêcT（uw 是常数，没有列）。第二遍是 ṙ=Buu̇+BgĠ：Bu 收集作用在两侧状态上的系数（cL+νp、cR−νp），Bg 收集作用在梯度上的系数（内部面在第 L、R 行写 ∓12νkT 到两个块，壁面面在第 L 行写 −νkT 到 GL）。没有替换，所以
As on the node-centred page, write the forward mode as a product of linear maps and transpose them one by one. Pass 1 is a 2N×N operator Γ (N=24 unknowns): Ġ=Γu̇, whose c-th block row is Σswcs(ês−êc)T, where the wall-midpoint term keeps only −wcsêcT (uw is a constant and has no column). Pass 2 is ṙ=Buu̇+BgĠ: Bu collects the coefficients on the two states (cL+νp, cR−νp), Bg those on the gradients (an interior face writes ∓12νkT into rows L and R in both blocks, a wall face −νkT into row L in the block GL). There is no replacement, so
用反向模式的记号：输入 w 就是残差的伴随 r̄，输出是状态的伴随 ū=ATr̄（代码里的 w 与 y），ū 与 Ḡ 先清零。从右往左读：先做第二遍的转置，对每个面取 h̄=r̄L−r̄R，把状态槽位直接累加进 ū，把梯度槽位存进 Ḡ=BgTr̄；再做第一遍的转置 ΓTḠ，把 Ḡ 经由权重送回 ū。壁面面在第一遍转置里只给 ūL 写 −wLR·ḠL：原来 uw 的那一份没有地方可去：ATw 里没有 uw 这一行，直接丢弃。第 11 节对壁值求导时，正是把这一份收集起来。
In reverse-mode notation the input w is the residual's adjoint r̄, and the output is the state's adjoint ū=ATr̄ (the code's w and y), with ū and Ḡ cleared first. Read from right to left: first the transpose of pass 2, taking h̄=r̄L−r̄R on every face, adding the state slots straight into ū and parking the gradient slots in Ḡ=BgTr̄; then the transpose of pass 1, ΓTḠ, sending Ḡ back into ū through the weights. In the transpose of pass 1 a wall face writes only −wLR·ḠL into ūL: the share that went to uw has nowhere to go, since ATw has no row for uw, and is dropped. The derivative with respect to the wall value in Section 11 is exactly this share, collected.
图 9.1 把三个函数逐项并排。与格点页面相比，这里没有“循环前清零、循环后加回”：弱施加不改动任何一行，转置也就不需要特别处理。
Figure 9.1 sets the three functions side by side. Compared with the node-centred page there is no "zero before the loop, add back after it": weak imposition changes no row, so the transpose needs no special treatment either.
| residual | residual_d | residual_b | |||||
|---|---|---|---|---|---|---|---|
| 执行顺序 | order | 第一遍 → 第二遍 | pass 1 → pass 2 | 第一遍 → 第二遍 | pass 1 → pass 2 | 第二遍T → 第一遍T | pass 2T → pass 1T |
| 第一遍 | pass 1 | Gc = Σs wcs(us − uc) | Ġc = Σs wcs(u̇s − u̇c) | ū += ΓTḠ （读 Ḡ，写 ū）(reads Ḡ, writes ū) | |||
| 第二遍 | pass 2 | r = Σ ± h | ṙ = Σ ± ḣ | h̄ = r̄L − r̄R → ū, Ḡ | |||
| 壁面 | the wall | uw 进入 hw 与 GL | uw enters hw and GL | uw 是常数，切线为 0 | uw is a constant: zero tangent | uw 的份额丢弃 | the share of uw is dropped |
图 9.2 逐句执行 residual_b。三行数组依次是输出 y（即 ū）、梯度的伴随 G_b（即 Ḡ）和输入 w（即 r̄）。第二遍倒着走时，每个面只把状态槽位写进 y，梯度槽位先存进 G_b；所有面走完之后，第一遍倒着走，才把 G_b 送回 y。
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̄). While pass 2 runs backwards each face writes only its state slots into y and parks its gradient slots in G_b; once every face is done, pass 1 runs backwards and sends G_b into y.
上面把两遍循环写成矩阵再转置。同一件事也可以在代码的层面做：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.
第一步，第二遍的最后两句（第 12、13 行）ṙL+=ḣ、ṙR−=ḣ，按规则 (c) 变成 residual_b 第 5 行：
Step 1, the last two statements of pass 2 (lines 12 and 13), ṙL+=ḣ and ṙR−=ḣ, become line 5 of residual_b by rule (c):
第二步，第 11 行 ḣ=(cL+νp)u̇L+(cR−νp)u̇R−νĠf·k，按规则 (b) 把 h̄ 分回三个输入；再把第 10 行 Ġf=12(ĠL+ĠR) 同样转置，就是 residual_b 第 6 至 8 行：
Step 2, line 11, ḣ=(cL+νp)u̇L+(cR−νp)u̇R−νĠf·k, hands h̄ back to its three inputs by rule (b); transposing line 10, Ġf=12(ĠL+ĠR), the same way gives lines 6 to 8 of residual_b:
状态直接累加进 ū（规则 (a)：原来读 u̇ 的地方现在往 ū 里累加）；梯度的份额先存进 Ḡ。壁面面上 Ġf=ĠL，系数 12 换成 1，也没有 ūR。
The states go straight into ū (rule (a): where u̇ was read, ū now accumulates); the gradients' share is parked in Ḡ. On a wall face Ġf=ĠL, the 12 becomes 1, and there is no ūR.
第三步，第一遍（第 4 至 6 行）Δ̇=u̇R−u̇L、ĠL+=wLRΔ̇、ĠR−=wRLΔ̇（Δ 是第 5 节第一遍跨过一个面的差值），倒着转置：先按 (b) 收集 Δ̄，再按 (a) 分给两个状态：
Step 3, pass 1 (lines 4 to 6), Δ̇=u̇R−u̇L, ĠL+=wLRΔ̇ and ĠR−=wRLΔ̇ (Δ 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 第 10 至 14 行的 γL−γR。
This is the γL−γR of lines 10 to 14 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Γ 里；壁面值 uw 是常数，在这两个矩阵里都没有列。
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Γ; the wall value uw is a constant, with no column in either matrix.
核对：用构建时生成的一对伪随机向量做点积测试，⟨w,Av⟩ 与 ⟨ATw,v⟩ 的相对差是 2.0×10−16；residual_b 与显式的 ATw 最多差 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 2.0×10−16 between ⟨w,Av⟩ and ⟨ATw,v⟩; residual_b differs from the explicit ATw 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.4）。
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 faces and design inputs included, in the order in which the reverse sweep actually runs (Figure 9.4).
第一步，一句的伴随。设一句算 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 第 5 至 8 行）。最后执行的分发最先转置：rL+=h 与 rR−=h 变成读 r̄；目标函数 J 是壁面面通量之和，所以带着 J̄ 扫描时，壁面面再加 J̄：
Step 2, pass 2 backwards, face by face (lines 5 to 8 of residual_b). The scatters ran last, so they are transposed first: rL+=h and rR−=h become reads of r̄; the objective J is the sum of the wall-face fluxes, so a sweep that carries J̄ adds J̄ on every wall face:
再转置通量这一句：两部分各按自己的偏导数把 h̄ 分回输入，几何量也在内：
Then the flux statement: each part hands h̄ back to its inputs with its own partial derivatives, the geometric ones included:
最后转置收集：平均梯度把 Ḡf 两侧各分一半，ḠL+=12Ḡf、ḠR+=12Ḡf（壁面面上 ḠL+=Ḡf）；右状态的伴随回到它的来处：内部面进右侧单元的 ū，壁面面进 ūw，入流面进 ū∞，出流面进 ūL（uR=uL，所以出流面上 ūL 共得 (cL+cR)h̄）。residual_b 只要 ū，把 ūw、ū∞、n̄、ξ̄ 这几份丢掉；第 11、12 节求设计导数时正是把它们收集起来。
Last, the gathers transposed: the averaged gradient gives each side half, ḠL+=12Ḡf and ḠR+=12Ḡf (ḠL+=Ḡf on a wall face); the right state's adjoint goes back where the state came from: into the right cell's ū on an interior face, into ūw on a wall face, into ū∞ on an inflow face, into ūL on an outflow face (there uR=uL, so ūL receives (cL+cR)h̄ in all). residual_b wants ū only, and drops the shares of ūw, ū∞, n̄ and ξ̄; the design derivatives of Sections 11 and 12 collect exactly those.
第三步，第一遍倒过来，逐面（第 10 至 14 行）。两句累加 GL+=wLRΔ、GR+=−wRLΔ 先收集 Δ̄，差值 Δ=uR−uL 再把它分给两个状态；权重的伴随留给第 12 节：
Step 3, pass 1 backwards, face by face (lines 10 to 14). The two accumulations GL+=wLRΔ and GR+=−wRLΔ first collect Δ̄, and the difference Δ=uR−uL then hands it to the two states; the weights' adjoints are for Section 12:
壁面面上 uR=uw，没有 GR：ūR 那一份进 ūw，residual_b 丢掉它，第 11 节对壁值求导时收集它。
On a wall face uR=uw and there is no GR: the share of ūR goes into ūw, which residual_b drops and Section 11 collects for the wall-value derivative.
第四步，源项。Rc←Rc−sc 转置成 s̄c−=r̄c；种子取 r̄=ψ 时就是第 11 节的 dJ/dsc=−ψc。
Step 4, the source. Rc←Rc−sc transposes into s̄c−=r̄c; seeded with r̄=ψ this is Section 11's dJ/dsc=−ψc.
图 9.4 把三列并排：左列是 residual 的语句，中列是切线（从上往下执行），右列是伴随（从下往上执行：右列自下而上读，就是 residual_b 自上而下的伪代码）。图 9.5 取面 f17（第 3 节画过的那个面，这张网格上最歪的内部面），把这些语句在本页 u* 上的数值全部写出来。第 8 节末尾（“本页一个面的算例：前向扫描”，图 8.4）与本节末尾（“同一个面的反向扫描”，图 9.9）两处推导取的是同一个面，种子相同，都是 v=du/du∞ 与 ψ，所以两处推导与这张表共有的数值一致；这张表多出几何量与通量的两个部分，两处推导多出整个模板的状态与逐级的内积。图 9.6 把 residual_b 整个跑一遍，种子取第 11 节的伴随解 ψ：ψ 满足 ATψ=−g，所以跑完以后 ū=−g，只在 16 个单元上不为零（8 个壁面单元和它们模板里的单元）。
Figure 9.4 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.5 takes the face f17 (the face drawn in Section 3, the most skewed interior face 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.9) work this same face, with the same seeds, v=du/du∞ and ψ, so the numbers they share with this table agree; this table adds the geometry and the two parts of the flux, and they add the states of the whole stencil and the inner products stage by stage. Figure 9.6 runs residual_b once through, seeded with the adjoint solution ψ of Section 11: ψ satisfies ATψ=−g, so at the end ū=−g, non-zero on 16 cells only (the 8 wall cells and the cells in their stencils).
核对：第一至四步逐句写成程序以后，用状态、来流值、源项、壁值、坐标上的随机种子与随机的 r̄、J̄，对偶恒等式 r̄Tṙ+J̄J̇=ūTu̇+D̄TḊ 两边的相对差最多 2.4×10−16；只对状态时与 residual_b 的相对差最多 1.9×10−16；种子取 r̄=ψ、J̄=1，一次扫描得到的 ū=ATψ+g 最大 4.6×10−16，D̄ 与本页代码分别算出的 dJ/du∞、dJ/ds、dJ/duw、dJ/dX 最多差 1.1×10−16；面 f17 上的对偶恒等式两边差 3.5×10−18；动画最后一帧与 −g 最多差 4.4×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 2.4×10−16 relative; on the state alone the program differs from residual_b by at most 1.9×10−16 relative; seeded with r̄=ψ and J̄=1, one sweep gives ū=ATψ+g no larger than 4.6×10−16, and a D̄ within 1.1×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 f17 differ by 3.5×10−18; and the animation's last frame differs from −g by at most 4.4×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 |
|---|---|---|---|---|---|---|---|
| 第一遍：梯度，逐个有 R 的面与壁面面 | pass 1: the gradients, every face with an R and every wall face | ||||||
| 4 | Δ=uR−uL 壁面面：uR=uw，ūw 丢弃wall face: uR=uw, ūw dropped | Δ̇=u̇R−u̇L | ūR+=Δ̄ ūL−=Δ̄ （第 11、12、14 行）(lines 11, 12, 14) | ||||
| 5 | GL+=wLRΔ | ĠL+=wLRΔ̇+ẇLRΔ | Δ̄+=wLR·ḠL w̄LR+=ΔḠL （第 10 行）(line 10) | ||||
| 6 | GR+=−wRLΔ | ĠR+=−wRLΔ̇−ẇRLΔ | Δ̄−=wRL·ḠR w̄RL−=ΔḠR （第 13 行）(line 13) | ||||
| 第二遍：通量，逐面 | pass 2: the fluxes, face by face | ||||||
| 9 | uL=u[L] uR=u[R] 边界面：uR 取 uw、u∞ 或 uLboundary face: uR is uw, u∞ or uL | u̇L=u̇[L] u̇R=u̇[R] | ū[L]+=ūL ū[R]+=ūR （第 6、7 行）(lines 6, 7) | ||||
| 10 | Gf=12(GL+GR) 壁面面：Gf=GLwall face: Gf=GL | Ġf=12(ĠL+ĠR) | ḠL+=12Ḡf ḠR+=12Ḡf （第 8 行）(line 8) | ||||
| 11 | ha=14vn(uL2+uR2)−ε|n|(uR−uL) 壁面面：ha=12uw2vn，没有耗散项，cL=0、cR=uwvnwall face: ha=12uw2vn, no dissipation term, cL=0 and cR=uwvn | ḣa=cLu̇L+cRu̇R+∂ha∂n·ṅ | ūL+=cLh̄ ūR+=cRh̄ n̄+=∂ha∂nh̄ （第 6、7 行）(lines 6, 7) | ||||
| 11 | hv=−ν[p(uR−uL)+Gf·k] 远场面：hv=0far field: hv=0 | ḣv=−ν[p(u̇R−u̇L)+Ġf·k]+∂hv∂n·ṅ+∂hv∂ξ·ξ̇ | ūL+=νph̄ ūR−=νph̄ Ḡf=−νkh̄ n̄,ξ̄+=… （第 6、7、8 行）(lines 6, 7, 8) | ||||
| 11 | h=ha+hv | ḣ=ḣa+ḣv | h̄a=h̄ h̄v=h̄ | ||||
| 12 | rL+=h | ṙL+=ḣ | h̄+=r̄L （第 5 行）(line 5) | ||||
| 13 | rR−=h | ṙR−=ḣ | h̄−=r̄R （第 5 行）(line 5) | ||||
| 循环之后：源项 | after the loops: the source | ||||||
| rc−=sc | ṙc−=ṡc | s̄c−=r̄c | |||||
| 量 | quantity | 式子 | formula | 数值 | value |
|---|---|---|---|---|---|
| primal：本页的 u* | primal: this page's u* | ||||
| uL,uR | 0.2613, 0.2506 | ||||
| n | (yb−ya,−(xb−xa)) | (−1.366, 0.04630) | |||
| |n|,vn | vn=β·n | 1.366, −1.342 | |||
| ξ | xR−xL | (−0.3195, 0.2882) | |||
| p | ξ·n/|ξ|2 | 2.429 | |||
| k | n−pξ | (−0.5895, −0.6536) | |||
| GL,GR | (−0.01891, −0.05842), (0.06801, 0.1460) | ||||
| Gf | 12(GL+GR) | (0.02455, 0.04381) | |||
| cL,cR | 12vnuL+ε|n|,12vnuR−ε|n| | 0.9176, −1.261 | |||
| ha,hv,h | −0.03219, 0.03466, 0.002473 | ||||
| 切线：种子是第 10 节的 v=du/du∞ | tangent: the seed is Section 10's v=du/du∞ | ||||
| u̇L,u̇R | 0.3573, 0.3501 | ||||
| ĠL,ĠR | (−0.03217, −0.06079), (0.08085, 0.2311) | ||||
| Ġf | 12(ĠL+ĠR) | (0.02434, 0.08513) | |||
| ḣa | cLu̇L+cRu̇R | −0.1136 | |||
| ḣv | −ν[p(u̇R−u̇L)+Ġf·k] | 0.04379 | |||
| ḣ | ḣa+ḣv | −0.06985 | |||
| 伴随：种子是第 11 节的 ψ | adjoint: the seed is Section 11's ψ | ||||
| r̄L,r̄R | ψL,ψR | −0.5340, −0.6008 | |||
| h̄ | r̄L−r̄R | 0.06678 | |||
| ūL | +=(cL+νp)h̄ | 0.1424 | |||
| ūR | +=(cR−νp)h̄ | −0.1653 | |||
| ḠL,ḠR | +=−12νkh̄ | (0.009842, 0.01091) | |||
| h̄ḣ | −0.004665 | ||||
| Σx̄ẋ | ūLu̇L+ūRu̇R+ḠL·ĠL+ḠR·ĠR | −0.004665 | |||
| 几何：第 12 节的反向扫描，这个面的种子 h̄f | geometry: Section 12's reverse sweep, this face's seed h̄f | ||||
| h̄f | ψL−ψR | 0.06678 | |||
| ∂h∂n | ∂ha∂n+∂hv∂n | (−0.001577, 0.006894) | |||
| ∂hv∂ξ | (0.03766, −0.003731) | ||||
| n̄ | h̄f∂h/∂n | (−1.05×10⁻⁴, 4.60×10⁻⁴) | |||
| ξ̄ | h̄f∂hv/∂ξ | (0.002515, −2.49×10⁻⁴) | |||
对任意 v、w，⟨w,Av⟩=⟨ATw,v⟩ 必须在舍入误差内成立。用第 8 节的记号，取 u̇=v、r̄=w，于是 ṙ=Au̇、ū=ATr̄，这就是反向模式的对偶恒等式 r̄Tṙ=ūTu̇。图 9.7 用两组随机向量检验。
For any v and w, ⟨w,Av⟩=⟨ATw,v⟩ must hold to round-off. In the notation of Section 8, with u̇=v and r̄=w, so that ṙ=Au̇ and ū=ATr̄, this is the duality identity of reverse mode, r̄Tṙ=ūTu̇. Figure 9.7 tests it on two random vectors.
第一步，测试的是什么。设 residual_b 算的是 Bw，B 是某个矩阵，写对了就是 AT。点积测试比较的两边之差是一个双线性型：
Step 1, what the test measures. Say residual_b computes Bw for some matrix B, which is AT if the code is right. The difference of the two sides the test compares is a bilinear form:
第二步，写对与写错。B=AT 时差只剩舍入误差。A−BT≠0 时，对固定的 w（只要 (A−BT)Tw≠0），wT(A−BT)v=0 只在一个超平面上成立，随机取的 v 落在上面的概率是零；差的大小与 A−BT 相对 A 的大小同一个量级。所以一对随机向量就足以判断：不用形成矩阵，也没有步长。
Step 2, right and wrong. With B=AT only round-off is left. With A−BT≠0, for a fixed w (as long as (A−BT)Tw≠0) the form wT(A−BT)v vanishes only on a hyperplane, which a random v hits with probability zero; its size is that of A−BT relative to A. So one random pair decides it: no matrix is formed, and there is no step size.
第三步，几种常见的写错，每种只改 residual_b 的一处（图 9.8）：
Step 3, some common slips, each changing one place of residual_b (Figure 9.8):
核对：每一列 30 对随机向量，residual_b 的各列用的是第 12 节表中那一组。本页的 residual_b 相对差最多 1.4×10−15（就是那张表里的数），几何的一对最多 6.5×10−15；几种写错里最小的也有 1.4×10−4，比舍入误差高出 10 个数量级以上（本页代码在构建时计算）。
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 1.4×10−15 relative (the number in that table), and the geometric pair within 6.5×10−15; the smallest error among the slips is 1.4×10−4, more than 10 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:
内部面通量向左加、向右减，所以先把这两条路径合并成通量的伴随 ；之后按通量、重构、原始状态的相反依赖顺序传播。反向语句中的 += 很重要：同一状态会被多个面和多个求积点读取，这些贡献必须累加。这里没有分部积分，也没有另行猜测连续伴随边界条件。
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 +=:
本页没有被替换的残差行：边界条件都在边界面的通量里，所以上面的面循环转置就是完整的 。
This page replaces no residual rows: the boundary conditions all sit in boundary-face fluxes, so the transposed face loop above is the complete .
一个面的切向量是四个初等 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:
仍从右往左读，这就是反向扫描的语句顺序：程序最后执行的分发最先处理，最先执行的取值最后处理。先把所有伴随置零，再以 为种子逐条执行：
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:
面循环结束后，源项的一对也按同样的规则反过来：，，。本页代码里每一条语句的伴随语句（顺序与前向相反）也都在本节里逐句写出。
After the face loop the source is reversed by the same rule: , , . 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:
这就是 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, , so that and are this face's contributions 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.534 | −0.004664501947 | |
| −0.6008 | |||
| （4′）取值(4′) gather | 0.06678 | −0.004664501947 | |
| （3′）通量的转置(3′) flux, transposed | 0.1424 | −0.004664501947 | |
| −0.1653 | |||
| 0.009842 0.01091 | |||
| 0.009842 0.01091 | |||
| （2′）（1′）加回(2′) (1′) add back | −0.02953 | −0.004664501947 | |
| 0.1913 | |||
| −0.1628 | |||
| 8.858×10⁻⁴ | |||
| −0.02275 | |||
| 0 |
4 级内积相同，最大相对差 ：对偶恒等式在这些数上成立；图 9.9 按反向语句逐条播放这些数。这个面不直接依赖 ，所以 ； 的伴随来自它真正进入的那些行。
The 4 inner products agree, with a largest relative difference of : the duality identity holds in these numbers. Figure 9.9 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.534， −0.6008。
(0) Zero and seed: , ; −0.534, −0.6008.
（4′）分发的转置是一次取值： 0.06678。内积 −0.004664501947。
(4′) The transposed scatter is a gather: 0.06678. Inner product −0.004664501947.
（3′）通量的转置：， 0.1424， −0.1653， (0.009842, 0.01091)， (0.009842, 0.01091)。内积 −0.004664501947。
(3′) The transposed flux: , 0.1424, −0.1653, (0.009842, 0.01091), (0.009842, 0.01091). Inner product −0.004664501947.
（2′）重构的转置：，，这个面给出 0（）。
(2′) The transposed reconstruction: and ; this face gives 0 ().
（1′）取值的转置是带累加的分发：， −0.02953， 0.1913， −0.1628， 8.858×10⁻⁴， −0.02275。内积 −0.004664501947。
(1′) The transposed gather is a scatter with accumulation: , −0.02953, 0.1913, −0.1628, 8.858×10⁻⁴, −0.02275. Inner product −0.004664501947.
对偶：4 级内积相同（最大相对差 9×10⁻¹⁶）， 在这个面上逐级成立。
Duality: the 4 inner products agree (largest relative difference 9×10⁻¹⁶); holds stage by stage on this face.
第 8 节只算了乘积 Av。前向模式真正要解的是切线方程：选设计变量的一个分量 Dk（种子 Ḋ=ek），对收敛解满足的 R=0 求导，u 随 Dk 的变化率 v=du/dDk，也就是这个种子下状态的切线 u̇，满足
Section 8 only computed the product Av. 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 R=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=∂R/∂u∞ 就是这些面上通量对右槽位的导数，写进各自的单元，共 4 个非零分量。
Here take the far-field free-stream value, Dk=u∞. Only the far-field inflow faces see it, as the right state of their convective flux (the far field carries no diffusive flux), so f=∂R/∂u∞ is the derivative of those fluxes with respect to their right slot, written into their cells: 4 non-zero entries.
切线方程是线性的，用第 6 节那个 Jacobi 迭代、同一个 diag(A) 求解：
The tangent equation is linear, and it is solved by the Jacobi iteration of Section 6, with the same diag(A):
每一轮只调用一次与矩阵有关的函数：第 8 节的 residual_d，算出 Av。diag(A) 是第 6 节的 diagonal 在 u* 处的值。除此之外，矩阵既不形成，也不分解：这就是 matrix-free 的前向求解。
Each sweep calls just one matrix-related function: residual_d of Section 8, which computes Av. diag(A) 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.
在解附近，primal 迭代的误差每一轮乘以 I−diag(A)−1A，切线迭代的误差乘的正是同一个矩阵。所以 primal 一过它非线性的起始阶段，两条收敛曲线就以同一个速率下降，这个速率就是谱半径 0.915：实测 k=150 到 250 之间，两者每轮都缩小 0.9149 倍（本页迭代矩阵第二大的特征值是 −0.86，衰减得慢，所以测量窗口比格点页面的 k=60 到 100 靠后）。起点不同：primal 从 u≡1 出发，切线从 v=0 出发，初始残差就是 f。图 10.1 逐轮播放切线迭代，看切线解 v 如何从远场入流面铺进区域。切线迭代在 k=323 时降到判据 10−13 以下，从 k=400 起落在舍入误差平台上。
Near the solution the primal iteration multiplies its error every sweep by I−diag(A)−1A, 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.915: measured between k=150 and 250, both shrink by 0.9149 per sweep (this page's iteration matrix has a second eigenvalue of −0.86, which fades slowly, so the window sits later than the node-centred page's k=60 to 100). 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 faces. The tangent iteration drops below the 10−13 criterion at k=323 and is on its round-off floor from k=400 on.
入流面上对流通量的右状态就是 u∞，所以 f 只来自这一项，写进它左侧的单元：
On an inflow face the right state of the convective flux is u∞ itself, so f comes from that term alone, written into the cell on its left:
入流面上 vn<0，两项都是负的：f 在 4 个单元上为负，其余为零。解 Av=−f 得到的 v 处处为正（图 10.2）：来流值抬高，整个区域的 u 都跟着抬高。
On an inflow face vn<0, so both terms are negative: f is negative in 4 cells and zero elsewhere. Solving Av=−f gives a v that is positive everywhere (Figure 10.2): raise the free-stream value and u rises over the whole domain.
核对：把 u∞ 挪动 ±10−6 重算残差（状态不动），中心差分与上式的 f 最多差 1.0×10−10（本页代码在构建时计算）。
Check: moving u∞ by ±10−6 and re-evaluating the residual (state held fixed), central differences differ from the f above by at most 1.0×10−10 (computed by this page's code at build time).
J 不直接依赖 u∞，所以
J does not depend on u∞ directly, so
本算例是 3.23751265643403；在 u∞±10−6 各重解一次 primal 的中心差分与它的相对差是 3.5×10−10。一次切线求解只回答一个设计分量：要得到对第 12 节 48 个坐标的导数，前向那一路要解 48 次。下一节从另一头走，一次伴随求解就得到全部。
Here it comes to 3.23751265643403; a central difference that re-solves the primal at u∞±10−6 differs from it by 3.5×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(A) 求解。迭代矩阵相同，收敛速率也相同，只是起点的残差 fk 不同：到判据要 254 到 323 轮。图 10.3 每类取一个，画出 f 与解出的 v。
Step 3, all are solved by the same Jacobi iteration with the same diag(A). The iteration matrix is the same, and so is the rate; only the starting residual fk differs: 254 to 323 sweeps to the criterion. Figure 10.3 takes one of each kind and draws f and the solution v.
核对：81 个分量各解一次，gTvk+∂J/∂Dk 与第 12 节一次伴随扫描给出的 D̄k 最多差 1.3×10−15（最大的分量是 3.238；本页代码在构建时计算）。
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 1.3×10−15 (the largest component is 3.238; 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 which only the paths and the source are active, and every of the iteration is one forward sweep with . All the work of a tangent solve is this face loop run again and again with different seeds, and the derivatives with respect to all design components take such solves.
伴随方程 ATψ=−g 的右端项 g=(∂J/∂u)T 由目标函数决定（第 7 节）。它同样用第 6 节的 Jacobi 迭代、同一个 diag(A) 求解（转置不改变对角元）：
The right-hand side of the adjoint equation ATψ=−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(A) (transposition leaves the diagonal alone):
与矩阵有关的唯一操作是一次第 9 节的 residual_b 调用：它的种子是 r̄=ψ，收敛时返回 ū=ATψ=−g。
The only matrix-related operation is one call to residual_b of Section 9: its seed is r̄=ψ, and at convergence it returns ū=ATψ=−g.
伴随解同样给出第 10 节的那个导数，不用再解任何方程：
The adjoint solution gives the derivative of Section 10 as well, with no further solve:
本算例两者都是 3.23751265643403，相对差 4.1×10−16。一次伴随求解回答一个输出对所有设计分量的导数：同一个 ψ 也给出第 12 节对 48 个坐标的导数，那里前向那一路要解 48 次。
Here both come to 3.23751265643403, with a relative difference of 4.1×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(A)−1AT，它是 I−Adiag(A)−1 的转置，而
The adjoint iteration multiplies its error every sweep by I−diag(A)−1AT, the transpose of I−Adiag(A)−1, and
与 primal、切线的 diag(A)−1A 相似。转置不改变特征值，所以三个迭代矩阵的谱完全相同：进入渐近阶段以后，primal、切线、伴随三条收敛曲线以同一个速率下降，这个速率就是谱半径 0.915。实测 k=150 到 250 之间，三条曲线每轮都缩小 0.9149 倍，四位有效数字相同。三条曲线的起点不同：primal 从 u≡1 出发，切线和伴随从 0 出发，初始残差就是 f 与 g；也只有 primal 有一段非线性的起始阶段。但一进入渐近阶段，它们就平行。图 11.1 把三条曲线画在一起，并逐轮播放伴随迭代，看伴随解 ψ 如何从壁面铺开：有扩散，上游、下游都有。
is similar to the primal's and the tangent's diag(A)−1A. 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.915. Measured between k=150 and 250, all three shrink by 0.9149 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=−Av，伴随方程给出 ATψ=−g；把前者代入、再用转置的定义挪一下括号：
Step 1, why ψTf=gTv. The tangent equation gives f=−Av and the adjoint equation ATψ=−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* 满足 Av*+f=0，所以 v*=v*−diag(A)−1(Av*+f)；与迭代式相减：
Step 2, the error of the tangent iteration. The exact solution v* satisfies Av*+f=0, so v*=v*−diag(A)−1(Av*+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(A) 是对角阵，(diag(A)−1)T=diag(A)−1，于是
Step 4, this matrix is similar to the transpose of the former. diag(A) is diagonal, (diag(A)−1)T=diag(A)−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 的相对差是 4.1×10−16；第四步的恒等式两边逐元素最多差 1.1×10−16，两个矩阵的 24 个特征值逐一配对最多差 3.6×10−15；第 150 到 250 轮实测，primal、前向、伴随每轮分别乘 0.91489、0.91489、0.91489，谱半径 ρ=0.91489（本页代码在构建时计算）。图 11.2 把三条曲线画到舍入误差平台。
Check: on this page ψTf and gTv differ by 4.1×10−16 relative; the two sides of the identity of step 4 differ entry by entry by at most 1.1×10−16, and the 24 eigenvalues of the two matrices, paired one by one, by at most 3.6×10−15; measured over sweeps 150 to 250, the primal, the tangent and the adjoint shrink per sweep by 0.91489, 0.91489 and 0.91489, against the spectral radius ρ=0.91489 (computed by this page's code at build time). Figure 11.2 draws the three curves down to their round-off floors.
收敛到机器精度是什么样子。三个迭代都没有停在判据处，而是跑满 600 轮。primal、切线、伴随分别在 k=315、323、321 时降到判据 10−13 以下，但这时仍以每轮 0.915 倍的比例下降：判据只是约定的停止点，并不说明迭代已经走到头。它们按同一比例一直降到 10−16 附近，然后分别从 k=386、400、400 起落在舍入误差平台上：此后再也不超过这次运行最后三分之一里的最大值，三个平台都在 1.6×10−16 与 5.5×10−16 之间。曲线变平，才说明迭代已经收敛到机器精度。停在判据处的 ψ 与平台上的 ψ 相差 3.7×10−13。
What convergence to machine precision looks like. None of the three iterations stops at the criterion; each runs the full 600 sweeps. The primal, the tangent and the adjoint drop below the 10−13 criterion at k=315, 323 and 321, but there they are still falling by a factor of 0.915 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=386, 400 and 400 on: from there they never again rise above the largest value of the run's last third, and the three floors all lie between 1.6×10−16 and 5.5×10−16. Only this flat tail shows that an iteration has converged to machine precision. The ψ at the criterion is 3.7×10−13 away from the ψ on the floor.
在单元 c 里加一个源项 sc，则 dJ/dsc=−ψc：−ψc 就是这个单位源最终流进物面的比例，其余的从远场流走（与中心差分的最大差为 4.4e−10）。壁面单元上是 0.70 到 0.74，内层三角形 0.53 到 0.62，外层三角形 0.41 到 0.60（图 11.3）。格点页面要把 ψ 在壁面节点上换成 −1 才得到这个场；这里物面上没有未知量，ψ 本身就是它。格点页面壁面上的那个 −1，在这里藏在目标函数本身：每个壁面面的种子是 1（第 12 节）。
Add a source sc in cell c; then dJ/dsc=−ψc, so −ψc is the fraction of that unit source that ends up in the body, the rest leaving through the far field (it agrees with central differences to 4.4e−10). It is 0.70 to 0.74 in the wall cells, 0.53 to 0.62 in the inner triangles and 0.41 to 0.60 in the outer ones (Figure 11.3). The node-centred page had to replace ψ by −1 on the wall nodes to get this field; here there are no unknowns on the body, and ψ is the field itself. The node page's −1 on the wall lives here in the objective: every wall face gets the seed 1 (Section 12).
第一步，强度为 sc 的源项进入单元 c 的平衡：Rc←Rc−sc（第 15 节的源项，一点求积）。它是设计变量 D 的一个分量，∂R/∂sc=−êc，而 J 不直接含它。于是由伴随公式：
Step 1, a source of strength sc enters the balance of cell c as Rc←Rc−sc (the source term of Section 15, one-point quadrature). It is one component of the design variables D, ∂R/∂sc=−êc, and J does not contain it directly. The adjoint formula gives
第二步，为什么是“份额”。把所有单元的平衡加起来：内部面的通量两两抵消，剩下物面、远场与源项。记 J∞ 为从远场流出的净通量，则在收敛解上
Step 2, why it is a share. Add up the balances of all cells: the fluxes through interior 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
一个单位源分成两份：−ψc 进物面，1+ψc 从远场出去。
A unit source splits in two: −ψc into the body and 1+ψc out through the far field.
核对：在 c1、 c10、 c9 各放一个 ±10−6 的源项并重解原始方程，dJ/ds 依次是 0.7442、0.6008、0.5340，与 −ψ 最多差 9.7×10−10；从远场出去的份额依次是 0.2558、0.3992、0.4660，两份之和与 1 最多差 1.6×10−9（本页代码在构建时计算）。
Check: a source of ±10−6 in c1, c10 and c9 in turn, re-solving the primal each time, gives dJ/ds= 0.7442, 0.6008 and 0.5340, within 9.7×10−10 of −ψ; the shares that leave through the far field are 0.2558, 0.3992 and 0.4660, and the two shares add up to 1 within 1.6×10−9 (computed by this page's code at build time).
壁值 uw 不是未知量，对状态 u 求导时它是常数（第 8 节）。但把每个壁面面上的壁值 uw,f 当作设计变量，照样能求 dJ/duw,f。它从两处进入：一是壁面面通量 hw，作为右状态；二是壁面单元 L 的最小二乘梯度，作为模板里壁面中点那一点，权重是 wLf（即第 9 节壁面面上的 wLR）。用第 12 节同一组面通量伴随 h̄f'=ψL−ψR+[f'∈Fw]，一次反向扫描就得到全部 8 个：
The wall value uw is not an unknown: differentiating with respect to the state u, it is a constant (Section 8). But take the wall value uw,f of each wall face as a design variable and dJ/duw,f follows all the same. It enters in two places: the wall-face flux hw, as its right state; and the least-squares gradient of the wall cell L, as the wall-midpoint point of the stencil, with weight wLf (the wLR of a wall face in Section 9). With the face-flux adjoints of Section 12, h̄f'=ψL−ψR+[f'∈Fw], one reverse sweep gives all 8:
第一项是壁面面通量 hw=12uw2vn−ν[p(uw−uL)+GL·k] 对右状态的直接导数；第二项经由最小二乘梯度，θf' 在内部面上取 12、在壁面面上取 1。第二项正是 residual_b（第 9 节）在第一遍转置里丢掉的那一份：它本该写进壁面中点这个“槽位”，而这个槽位不是未知量；要对壁值求导，就把它收集起来。本算例的 8 个值介于 −0.312 与 −0.126 之间：提高某个壁面面的壁值，流进物面的通量就减少，与格点页面的 dJ/duw,i=−ψi 同号。经由最小二乘梯度的部分只占 0.6% 到 3.9%，其余来自两点部分 −νph̄f（uw=0，对流部分的导数 uwvn 为零）。与中心差分（每次重解原始方程）的最大差是 3.5e−10。所以弱施加并不妨碍对壁值求导，只是导数不像格点页面那样直接写在 ψ 的某个分量上。
The first term is the direct derivative of the wall-face flux hw=12uw2vn−ν[p(uw−uL)+GL·k] with respect to its right state; the second goes through the least-squares gradient, with θf' equal to 12 on an interior face and 1 on the wall face. The second term is exactly the share that residual_b (Section 9) drops in its transposed pass 1: it belongs to the wall-midpoint slot, which is not an unknown, and to differentiate with respect to the wall value one collects it instead. Here the 8 values lie between −0.312 and −0.126: raising the wall value on one wall face lowers the flux into the body, the same sign as the node-centred page's dJ/duw,i=−ψi. The part through the least-squares gradient is only 0.6% to 3.9%; the rest is the two-point part −νph̄f (with uw=0 the convective derivative uwvn vanishes). The largest difference from central differences (re-solving the primal each time) is 3.5e−10. So weak imposition does not stand in the way of wall-value derivatives; they are just not written directly in one component of ψ, as on the node-centred page.
−ψc 是单元 c 里一个单位源流进物面的份额（上文）。源就放在物面上时，它全部被吸收，份额是 1：这正是附录 A.1 第 3 小节连续伴随的壁面条件 λ=1（附录 A.1 第 1 小节定义 λh=−ψh，附录 A.3 用它与离散伴随比较）。本页粗网格上，贴壁单元的中心离物面还有一段距离，−ψ 只有 0.70 到 0.74；加密网格，贴壁的值向 1 靠拢。
−ψc is the share of a unit source in cell c that ends up in the body (above). A source placed on the body is absorbed entirely, a share of 1: 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). On this page's coarse mesh the centres of the wall cells sit some way from the body, and −ψ is only 0.70 to 0.74 there; refine the mesh and the values next to the wall close in on 1.
伴随逆着来流输运。伴随方程的矩阵是 AT。在 A 里，一个面的对流通量主要读迎风一侧的状态，一个状态的变化顺着 β 传到下游的行；转置把这种依赖倒过来，伴随在一点上的值听的是它下游的伴随，于是伴随逆着 β 传。连续问题里也是这样：线性化的对流 uβ·∇(·) 转置以后是 −∇·(uβ·)，方向相反。而 u 在物面附近很小（物面上是 0），所以贴壁的一层里伴随主要靠扩散，离开物面、u 变大以后，逆流的输运才显出来。primal 正好相反：u 顺着 β 从远场带进来，迎流一侧在物面前形成一层薄的边界层，背流一侧留下长长的尾迹。
The adjoint is transported against the flow. The matrix of the adjoint equation is AT. In A 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.719 对 0.601、0.781 对 0.674、0.808 对 0.714、0.835 对 0.727：每一层上游都更大，即伴随逆着来流伸得更远。u 则是每一层物面的影响都顺着来流伸得更远：第 3 层上游已恢复到 0.659，下游只有 0.120。沿上游射线，贴壁的值在第 0 至 3 层依次是 0.719、0.872、0.953、0.981，向 1 靠拢。
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.719 against 0.601, 0.781 against 0.674, 0.808 against 0.714 and 0.835 against 0.727: the upstream value is larger on every level, so the adjoint reaches further against the flow. For u the body's influence reaches further along the flow on every level: on level 3 it has recovered to 0.659 upstream and is only 0.120 downstream. Along the upstream ray the value next to the wall is 0.719, 0.872, 0.953, 0.981 on levels 0 to 3, closing in on 1.
核对：每一层存下的 ψ 代回伴随方程，ATψ+g 的均方根最多 1.9×10−16；ψTf 与同一层切线解的 gTv（附录 A.3 表中的 dJ/du∞）相对差最多 4.1×10−16（本页代码在构建时计算）。
Check: put back into the adjoint equation, each level's stored ψ leaves ATψ+g with an rms of at most 1.9×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.1×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 | 次，系数 , 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 阶的线性方程组，伴随路线只解 1 个。设计变量多而目标少时（形状优化通常如此）伴随路线便宜得多；设计变量少而目标多时，前向路线更合适。
This page has and : the tangent route solves 81 linear systems of order 24, 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∂R/∂X。用反向模式的记号，这就是坐标的伴随 X̄=(∂R/∂X)Tψ+(∂J/∂X)T，种子是残差的伴随 r̄=ψ 与目标函数的伴随 J̄=1。ψ 那一项是残差在状态固定时对坐标的反向扫描，每个面的种子是 ψL−ψR；J 那一项来自 J=Σf∈Fwhw 本身对坐标的依赖，也就是每个壁面面再加一个种子 1。两项合在同一次扫描里，面通量的伴随 h̄f 就是两者之和：
dJ/dX=∂J/∂X+ψT∂R/∂X. In reverse-mode notation this is the adjoint of the coordinates, X̄=(∂R/∂X)Tψ+(∂J/∂X)T, seeded with the residual's adjoint r̄=ψ and the objective's J̄=1. The ψ term is a reverse sweep of the residual over the coordinates with the state held fixed, each face seeded with ψL−ψR; the J term comes from J=Σf∈Fwhw depending on the coordinates itself, which is one more seed of 1 on every wall face. Both go into the same sweep, and the adjoint h̄f of a face flux is their sum:
一个面的通量通过这些几何量依赖坐标（图 12.1）：面法向 n=(yb−ya,−(xb−xa))，由面的两个端点 a、b 直接给出；各单元的顶点平均点 xc 与壁面面中点，由它们得到 ξ，再得到 p、k；以及经由模板向量进入的最小二乘权重。下面逐条写出每条路径上的偏导，再把它们串成一次反向扫描。
A face's flux depends on the coordinates through these geometric quantities (Figure 12.1): the face normal n=(yb−ya,−(xb−xa)), given directly by the face's two end nodes a and b; the cells' vertex averages xc and the wall-face midpoints, which give ξ, and from ξ the p and k; and the least-squares weights, through the stencil vectors. Below, the partial derivatives along each path are written out one by one and then strung together into one reverse sweep.
扩散通量对 n 与 ξ。用 n=pξ+k，即 Gf·k=Gf·n−pGf·ξ，把扩散通量写成只含 n 与 ξ 的形式：
The diffusive flux, with respect to n and ξ. With n=pξ+k, that is Gf·k=Gf·n−pGf·ξ, the diffusive flux is written in n and ξ alone:
状态与梯度固定时，p 对 n 的偏导是 ξ/|ξ|2，对 ξ 的偏导是 (n−2pξ)/|ξ|2，于是
With the states and the gradients held fixed, p has the partial derivative ξ/|ξ|2 with respect to n and (n−2pξ)/|ξ|2 with respect to ξ, so that
壁面面上 Δu=uw−uL，Gf=GL。对流部分只经由 n：内部面与远场面上 ∂ha/∂n=14(uL2+uR2)β−ε(uR−uL)n/|n|（远场的入流、出流划分保持冻结）；壁面面上 ha=12uw2vn，∂ha/∂n=12uw2β，本页 uw=0，它为零。
On a wall face Δu=uw−uL and Gf=GL. The convective part depends on n only: on interior and far-field faces ∂ha/∂n=14(uL2+uR2)β−ε(uR−uL)n/|n| (with the far field's inflow/outflow split kept frozen); on a wall face ha=12uw2vn and ∂ha/∂n=12uw2β, zero here since uw=0.
最小二乘梯度对模板向量。单元梯度是 Gc=Mc−1bc，其中 bc=Σsξs(us−uc)，Mc=ΣsξsξsT；壁面中点那一点的 us 就是 uw。状态固定，对模板向量做前向求导，ξ̇s 是它们的切线：
The least-squares gradient, with respect to the stencil vectors. A cell gradient is Gc=Mc−1bc with bc=Σsξs(us−uc) and Mc=ΣsξsξsT, where us at a wall midpoint is uw. Holding the state fixed, differentiate in forward mode with respect to the stencil vectors, whose tangents are ξ̇s:
反向扫描要的是：梯度的伴随 Ḡc 怎样变成各模板向量的伴随 ξ̄s。它们由对偶恒等式 Ḡc·Ġc=Σsξ̄s·ξ̇s 定义：把上式代入左边，用 Mc 的对称性把 Mc−1 移到 Ḡc 上，令 b̄=Mc−1Ḡc（即 bc 的伴随），再按 ξ̇s 整理：
What the reverse sweep needs is how the gradient's adjoint Ḡc becomes the adjoints ξ̄s of the stencil vectors. They are defined by the duality identity Ḡc·Ġc=Σsξ̄s·ξ̇s: substitute the above into its left side, move Mc−1 onto Ḡc by the symmetry of Mc, set b̄=Mc−1Ḡc (the adjoint of bc), and collect the terms by ξ̇s:
这里的 Ḡc 由平均梯度那一路收集：每个面把 −θνkh̄f 加给它两侧的单元（壁面面只有一侧），θ 在内部面上是 12，在壁面面上是 1。这与第 9 节的 Ḡ 是同一个公式，只是种子换成了这次扫描的面通量伴随 h̄f。
The Ḡc here is collected along the averaged-gradient path: each face adds −θνkh̄f to the cells on its two sides (a wall face has one), with θ=12 on an interior face and 1 on a wall face. It is the formula of Section 9's Ḡ, with the face adjoints h̄f of this sweep as the seed.
串成一次反向扫描。每个面先得到 n̄=h̄f∂h/∂n 与 ξ̄=h̄f∂hv/∂ξ；再由各单元的 Ḡc 得到各模板向量的 ξ̄s。ξ 与 ξs 都是两点之差，终点的伴随加上它，起点的伴随减去它。最后分给节点：顶点平均点那一端平均分给单元的各个顶点（三角形每个 13，四边形每个 14），壁面面中点那一端各分一半给面的两个端点，n̄ 按 n=(yb−ya,−(xb−xa)) 直接交给面的两个端点。
Strung into one reverse sweep. Each face first gets n̄=h̄f∂h/∂n and ξ̄=h̄f∂hv/∂ξ; then each cell's Ḡc gives its stencil vectors their ξ̄s. Both ξ and ξs are differences of two points, so the end point's adjoint gains ξ̄ and the start point's loses it. Finally the points hand it to the nodes: a vertex-average end is shared equally among the cell's vertices (13 each on a triangle, 14 on a quadrilateral), a wall-midpoint end half-and-half between the face's two end nodes, and n̄ goes straight to the face's two end nodes through n=(yb−ya,−(xb−xa)).
一个坐标分量 Xm 并不牵动所有面。在 图 12.2 里点选一个网格节点，可以看到它移动时跟着变的几何量：法向 n 变化的面（紫），移动的顶点平均点与壁面中点（琥珀色点），ξ 从而 p、k 变化的面（虚线），最小二乘权重变化的单元（着色），以及通量因此变化的面（红）。反向扫描求 dJ/dXm 时，只有这些红色的面有贡献，每个面带着自己的面通量伴随 h̄f（标在面上）；紫色箭头就是这样得到的 dJ/dXm。
One coordinate Xm does not reach every face. Pick a mesh node in Figure 12.2 to see what moves with it: the faces whose normal n moves (violet), the vertex averages and wall midpoints that move (amber dots), the faces whose ξ, and so p and k, move (dashed), the cells 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̄f (written on it); the violet arrow is the dJ/dXm they give.
第一步，依赖什么。Gc=Σswcs(us−uc) 对状态是线性的；权重 wcs=Mc−1ξs 只依赖几何（不加距离权，与 u 无关），经模板向量 ξs 依赖坐标。所以对状态求切线时权重不动，对坐标求切线时才要权重的切线。
Step 1, what depends on what. Gc=Σswcs(us−uc) is linear in the state; the weights wcs=Mc−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.
第二步，权重的切线。由 McMc−1=I 得 ddtMc−1=−Mc−1ṀcMc−1，于是
Step 2, the weights' tangent. From McMc−1=I, ddtMc−1=−Mc−1ṀcMc−1, and so
第三步，梯度的切线有两项：
Step 3, the gradient's tangent has two terms:
第一项是第 8 节的第一遍；第二项只在坐标带切线时出现。其中 ḃc=Σsξ̇s(us−uc) 是 bc=Σsξs(us−uc) 在状态不动时的切线，只含几何的切线；代入第二步，第二项就是上文状态固定时的 Ġc。
The first term is pass 1 of Section 8; the second appears only for a coordinate seed. In it ḃc=Σsξ̇s(us−uc) is the tangent of bc=Σsξs(us−uc) with the state held fixed, the geometric tangents only; with step 2, the second term is the Ġc at a fixed state above.
第四步，对线性场精确，网格移动时也不丢。对 ΣswcsξsT=I 求切线，得 Σs(ẇcsξsT+wcsξ̇sT)=0。线性场 u=u0+∇u·x 在移动的点上取值时 us−uc=ξs·∇u、u̇s−u̇c=ξ̇s·∇u，代入第三步：
Step 4, exactness on linear fields survives a moving mesh. The tangent of ΣswcsξsT=I is Σs(ẇcsξsT+wcsξ̇sT)=0. A linear field u=u0+∇u·x sampled at moving points has us−uc=ξs·∇u and u̇s−u̇c=ξ̇s·∇u, and step 3 gives
网格怎么动，线性场的梯度都还是 ∇u。
However the mesh moves, a linear field's gradient stays ∇u.
第五步，转置：梯度的伴随经同一组权重回到模板的各点。把第三步配成对偶恒等式 Ḡc·Ġc=Σsūsu̇s+ūcu̇c+Σsw̄cs·ẇcs，按各个切线读出系数：
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 Ḡc·Ġc=Σsūsu̇s+ūcu̇c+Σsw̄cs·ẇcs and read off the coefficients of the tangents:
前两式是“转置的模板”：前向时 Gc 从模板各点收集，伴随时 Ḡc 按同样的权重分回这些点，每点一份，自己减去总和。这就是 residual_b 第 10 至 14 行的 γ。壁面中点那一份是 ūw：residual_b 丢掉它，第 11 节收集它。
The first two are the transposed stencil: forward, Gc gathers from the stencil's points; in the adjoint, Ḡc hands back to the same points with the same weights, one share each, and the cell itself gives up their sum. These are the γ of lines 10 to 14 of residual_b. The wall midpoint's share is ūw: residual_b drops it, Section 11 collects it.
第六步，再把第二步转置，从权重的伴随得到模板向量的伴随：
Step 6, transpose step 2, and the weights' adjoints give the stencil vectors':
代入第五步的 w̄cs=Ḡc(us−uc)，并记 b̄=Mc−1Ḡc：Σsw̄csξsT=ḠcbcT，所以 M̄c=−b̄GcT，三项合起来正是上文的 ξ̄s=b̄(us−uc)−b̄(Gc·ξs)−Gc(b̄·ξs)。
Put in step 5's w̄cs=Ḡc(us−uc) and write b̄=Mc−1Ḡc: then Σsw̄csξsT=ḠcbcT, so M̄c=−b̄GcT, and the three terms together are exactly the ξ̄s=b̄(us−uc)−b̄(Gc·ξs)−Gc(b̄·ξs) above.
第七步，一个单元的数值：壁面单元 c1（4 个模板点，其中一个是壁面中点），Ḡc 取本节这一次扫描中收集到的 (0.0191,−0.0173)，b̄=Mc−1Ḡc=(0.0507,−0.0186)：
Step 7, the numbers of one cell: the wall cell c1 (4 stencil points, one of them the wall midpoint), with the Ḡc=(0.0191,−0.0173) collected by this section's sweep and b̄=Mc−1Ḡc=(0.0507,−0.0186):
| 模板点 s | stencil point s | ξs | wcs | us−uc | wcs·Ḡc | ξ̄s |
|---|---|---|---|---|---|---|
| c2c2 | (−0.4866, 0.4269) | (−1.295, 0.4558) | 0.08300 | −0.03261 | (0.005699, 2.94×10⁻⁴) | |
| c8c8 | (−0.03617, −0.6634) | (−0.2128, −0.8411) | 0.06132 | 0.01053 | (0.004577, −0.002455) | |
| c10c10 | (0.2802, 0.4099) | (0.8560, 0.5629) | 0.07193 | 0.006562 | (0.001503, −0.001034) | |
| f33 的中点the midpoint of f33 | (−0.2055, −0.1025) | (−0.5945, −0.1634) | −0.1786 | −0.008505 | (−0.007867, 0.003517) |
核对：权重的切线与第 14 节的复数步长最多差 8.9×10−16；线性场放在随机移动的网格上，梯度的切线最大 4.4×10−16（第四步）；第五、六步包含在第 9 节的对偶恒等式核对里（本页代码在构建时计算）。
Check: the weights' tangent differs from the complex step of Section 14 by at most 8.9×10−16; for a linear field on a randomly moving mesh the gradient's tangent is at most 4.4×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̄f 经 ∂h/∂n、∂hv/∂ξ 变成 n 与 ξ 的伴随，推动面的端点和 ξ 的两端（顶点平均点再平均分给各顶点，壁面中点分给两个端点），梯度的份额存进 Ḡ。所有面走完、Ḡ 齐了以后，再逐单元：Ḡc 经最小二乘权重变成模板向量的伴随（上一段第六步），推动模板的各点。代码先把各点的伴随收集起来、最后再分给节点；动画每一步直接分，和是一样的。每一步的推动都成对出现、大小相等方向相反（n、ξ、ξs 都是两点之差），所以任何时刻全部箭头之和都停在舍入误差上：“整张网格平移不改变任何东西”逐面成立。最后一帧就是下面 图 12.4 的箭头。
Figure 12.3 plays this reverse sweep step by step. Face by face first: each face turns its adjoint h̄f, through ∂h/∂n and ∂hv/∂ξ, into adjoints of n and ξ that push the face's end nodes and the two ends of ξ (a vertex average shares its push equally among its vertices, a wall midpoint between the face's two end nodes), and parks its gradient share in Ḡ. Once every face is done and Ḡ is complete, cell by cell: Ḡc 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；每一步推动的和最大 2.8×10−17，全部箭头之和在任何一帧最大 1.9×10−16（本页代码在构建时计算）。
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 2.8×10−17, and the sum of all arrows in any frame to at most 1.9×10−16 (computed by this page's code at build time).
图 12.4 画出这一次反向扫描的结果：每个网格节点上一个箭头 dJ/dXi。物面节点上的 8 个箭头全部朝外（与径向最多差 7.3°）：把物面往外推，物体吸收的通量变大。
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 7.3° of the radial direction): push the body surface out, and the flux the body absorbs grows.
整张网格平移不改变任何东西，所以 ΣidJ/dXi 在两个方向上都必须为零：实测为 0 与 1.4e−16。
Translating the whole mesh changes nothing, so ΣidJ/dXi must vanish in both directions: measured, 0 and 1.4e−16.
| 精确 Jacobian 对中心差分（相对） | exact Jacobian against central differences (relative) | 4.3e−11 |
|---|---|---|
| residual_d、residual_b 对显式 A（第 9 节的那一对随机向量） | residual_d, residual_b against the explicit A (the random pair of Section 9) | 4.4e−16, 1.3e−15 |
| 点积测试（30 组随机向量中最差） | dot-product test (worst of 30 random pairs) | 1.4e−15 |
| Jacobi 满足判据的轮次（primal / 切线 / 伴随） | Jacobi sweeps to the criterion (primal / tangent / adjoint) | 315 / 323 / 321 |
| dJ/du∞：切线 gTv 对伴随 ψTf（相对） | dJ/du∞: tangent gTv against adjoint ψTf (relative) | 4.1e−16 |
| dJ/dX 对全链路中心差分（h = 10⁻⁶，每次重解） | dJ/dX against full-chain central differences (h = 10⁻⁶, re-solving each time) | 绝对 1.5e−9，相对 3.7e−9absolute 1.5e−9, relative 3.7e−9 |
| Σi dJ/dXi | Σi dJ/dXi | (0, 1.4e−16) |
| dJ/dsc = −ψc 对中心差分 | dJ/dsc = −ψc against central differences | 4.4e−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̄c=−ψc（第 11 节），壁值的 ūw（第 11 节），坐标的 X̄（本节）。
Its four parts all appear above: ū∞=ψTf (Section 11), s̄c=−ψc (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 次求解各自到判据，共 23431 轮（每轮一次 residual_d）；伴随：一次求解 321 轮（每轮一次 residual_b），再加一次反向扫描。求解次数之比是 81 比 1，轮数之比约 73 比 1。图 12.5 把两条路并排：伴随一步就立起全部 81 根柱，前向一个一个地把绿点落在同一根柱上。
Step 4, the cost. Everything that involves the matrix sits in a Jacobi iteration, one matrix product per sweep. Tangent: 81 solves, each to the criterion, 23431 sweeps in all, one residual_d each; adjoint: one solve of 321 sweeps, one residual_b each, and one reverse sweep. The solves are 81 against 1, the sweeps about 73 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.
核对：两组 81 个数最多差 1.3×10−15（最大的分量是 3.238；本页代码在构建时计算）。
Check: the two sets of 81 numbers differ by at most 1.3×10−15 (the largest component is 3.238; 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、最小二乘权重），重新解一次非线性的原始方程，再把 8 个壁面面的通量加起来。所以这里验证的是全导数，状态对几何的响应也算在内。代价也正因如此：中心差分每个分量要两次非线性求解，一共 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 (face normals, vertex averages, ξ, p, k, least-squares weights), solves the nonlinear primal again, and only then sums the fluxes through the 8 wall faces. 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 个分量逐一扫描：中心差分的最优步长有 33 个落在 h=10−5，其余 15 个在 10−4；前向差分有 37 个落在 10−7，其余在 10−6 或 10−8。这与 ϵmach1/3、ϵmach1/2 的估计（ϵmach 为机器精度）一致。各分量在最优步长下能达到的相对误差，中心差分的中位数是 1.1e−9，最好 2.5e−12；前向差分的中位数是 1.0e−7。默认的 h=10−6 对中心差分略偏小，对前向差分略偏大：在这个步长下，48 个分量中差分与伴随的最大绝对差，中心差分是 1.5e−9，前向差分是 4.2e−7。点积测试（第 9 节）则完全不含步长，是更干净的检查。差分的步长两难，在第 14 节的复数步长那里不存在。
Where the optimal step lies. Sweep the 48 components one by one: the central optimum is h=10−5 for 33 of them and 10−4 for the other 15; the forward optimum is 10−7 for 37 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 1.1e−9 for central differences, 2.5e−12 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.5e−9 for central differences and 4.2e−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 中最大的分量（节点 7 的 y 坐标），步长从 10−1 到 10−12 每次重解原始方程：中心差分在 h=10−5 处相对误差最小，为 2.3×10−10；前向差分在 h=10−7 处最小，为 2.8×10−8（本页代码在构建时计算）。
Check: Figure 13.2 takes the largest component of dJ/dX (the y coordinate of node 7) and re-solves the primal at every step from 10−1 to 10−12: central differences are best at h=10−5, with a relative error of 2.3×10−10; forward differences at h=10−7, with 2.8×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 矩阵的逆）都要重算一遍。本页用到复数的加、减、乘、除，以及一次开方：|n|=n·n，取主值平方根。远场面是入流还是出流，沿用实数程序在基准网格上的判断，不随虚部变。两遍循环与 Jacobi 迭代（diag(A) 也由同一个复数循环收集）都不变，只是每个数多带一个虚部；判据要连虚部一起看：ImR/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 vertex averages, the face vector n and |n|, ξ, p, k and the least-squares weights (with the inverse of their 2×2 matrix) are all recomputed. 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(A) 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 ImR/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* 处展开，R(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 R(u*)=0:
于是 Jacobi 更新 u←u−diag(A)−1R(u) 的虚部除以 h 是
so the imaginary part of the Jacobi update u←u−diag(A)−1R(u), divided by h, is
这里 diag(A) 也是复数（它由同一个复数循环收集），所以更新里那个商的虚部完整写出来是
Here diag(A) is complex as well (the same complex loop collects it), so the imaginary part of the quotient in the update reads in full
从 u* 出发，ReR 停在舍入误差上，第二项也只有舍入误差的量级，剩下的正是第 10 节的前向 Jacobi 迭代：v0=0，同一个 diag(A)，右端项 f=∂R/∂Xm 取自第 12 节的几何反向代码（逐行取出，是手写的导数，不是复数步长）。所以每一轮的 Imuk/h 就是前向迭代的 vk，不只是收敛后一致：48 个坐标分量、600 轮中的每一轮，最大差（除以 max|v|）是 5.6×10−15。实部与实数的 Jacobi 迭代相差不超过 1.1×10−16，但不是逐比特相同：实数程序用 Math.hypot 算 |n|，复数程序用平方根，两者在 48 个面中的 14 个上末位舍入不同。从初始流场出发就不同了：ReR≠0，第二项留下来，迭代把 diag(A) 也求了导，虚部是 primal 迭代值对坐标的导数：另一个序列，只在收敛时才到同一个极限。前几轮的差可达 max|v| 的 2.3 倍，而且每个分量都要把 primal 重解一遍：节点 1 的 x 从初始流场出发，虚部的残差要到第 333 轮才满足判据；从 u* 出发与前向迭代一样，是第 296 轮。
Started at u*, ReR 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(A), and the right-hand side f=∂R/∂Xm 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 600 sweeps the largest difference (over max|v|) is 5.6×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 14 of the 48 faces. Started at the initial field it is different: ReR≠0, the second term stays, the iteration differentiates diag(A) 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 2.3 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 333; started at u*, at sweep 296, the forward iteration's count.
每个坐标分量一次复数求解，得到梯度的一个分量。48 个分量与伴随梯度之差的最大值（除以梯度的最大分量）是 6.7×10−16，第 13 节的中心差分（h=10−6）是 3.7×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 6.7×10−16, against 3.7×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:
其中 |Ωc| 是单元 c 的面积，xc 是它顶点的平均点。
where |Ωc| is the area of cell c and xc the average of its vertices.
壁面：uw=ue，在壁面面的中点求值，照样弱施加：它同时进入壁面面的通量与壁面单元的最小二乘模板。远场：入流面的外侧状态取 ue，在面中点 xf 求值；原问题里为零的远场扩散通量，换成精确解的法向通量：
Wall: uw=ue at the wall-face midpoints, imposed weakly as before: it enters both the wall-face flux and the wall cells' least-squares stencils. Far field: the outer state on an inflow face is ue at the 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.
误差在单元上量，与 ue 在顶点平均点 xc 的值比较，用单元面积加权的 L2 范数与最大范数：
The error is measured on the cells, against ue at the vertex average xc, in the L2 norm weighted by the cell 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 越慢：只有扩散时，这一列最细的第 3 层要 13012 轮才满足判据；完整方程最细的第 5 层要 96264 轮。从第 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: with diffusion alone the finest level of its column, level 3, takes 13012 sweeps to meet the criterion, and the full equation's finest, level 5, takes 96264. 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 ℓ | 单元数 | cells | 只有扩散（β = 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 | 6.88e−3 | — | 1.38e−2 | — | 2.80e−2 | — | 5.42e−2 | — | ||||||||||
| 1 | 96 | 1.80e−3 | 1.94 | 5.34e−3 | 1.37 | 1.86e−2 | 0.59 | 3.65e−2 | 0.57 | ||||||||||
| 2 | 384 | 4.42e−4 | 2.02 | 1.40e−3 | 1.93 | 1.14e−2 | 0.71 | 2.37e−2 | 0.62 | ||||||||||
| 3 | 1536 | 1.10e−4 | 2.01 | 3.60e−4 | 1.96 | 6.42e−3 | 0.82 | 1.36e−2 | 0.80 | ||||||||||
| 4 | 6144 | Jacobi 迭代发散Jacobi diverges | 3.43e−3 | 0.90 | 7.33e−3 | 0.89 | |||||||||||||
| 5 | 24576 | Jacobi 迭代发散Jacobi diverges | 1.78e−3 | 0.95 | 3.81e−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.94、2.02、2.01，最大范数在第 3 层也到了 1.96；格点页面对应的 L2 阶是 1.99、2.09、2.08、2.05、2.03（算到第 5 层）。二阶的前提是相容：最小二乘梯度与带修正的平均梯度通量对线性场都精确（第 3 节）。测得的阶说明，单点积分的源项与不规则的网格都没有把它拉低。
With β=0 and ε=0 only diffusion and the source remain. The L2 orders on this page are 1.94, 2.02 and 2.01, and the max-norm order reaches 1.96 on level 3; the node-centred page's L2 orders are 1.99, 2.09, 2.08, 2.05 and 2.03 (to level 5). 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.
这一列为什么只到第 3 层。从第 4 层起，本页的 Jacobi 迭代在这个问题上发散。只有扩散时残差对 u 是仿射的，R(u)=Au−b，一轮迭代把误差乘以 M=I−diag(A)−1A；只要 M 有一个特征值落在单位圆外，迭代就发散。用幂迭代量出 M 最负的特征值，第 0 层到第 5 层依次是 −0.917、
Why this column stops at level 3. From level 4 on, the page's Jacobi iteration diverges on this problem. With diffusion alone the residual is affine in u, R(u)=Au−b, and one sweep multiplies the error by M=I−diag(A)−1A; as soon as one eigenvalue of M lies outside the unit circle, the iteration diverges. The most negative eigenvalue of M, measured by power iteration, is −0.917, −0.935, −0.973, −0.993, −1.004 and −1.007 on levels 0 to 5: it crosses −1 at level 4. The corresponding mode changes sign every sweep (the eigenvalue is negative), has opposite signs across 87.8% of the interior faces, and has almost all its energy in the two cell layers next to the far field (94.6% on level 4); in this problem the far-field faces carry a prescribed diffusive flux, which contributes nothing to the diagonal. On level 4 the iteration started from ue (every level of the study starts from ue sampled at the unknowns) eventually multiplies its residual by 1.0036 per sweep (the average over the last 5000 sweeps), the magnitude of that eigenvalue; in 10000 sweeps the residual grows from 2.23e−5 to 2.03e11. Every solve on this page is the same Jacobi iteration, so this column stops at level 3; the full equation (whose diagonal also holds the convective and dissipative shares) converges on the same family all the way to level 5.
对流通量里的耗散项可以写成扩散两点部分的样子：
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.487、0.270、0.140、0.071、0.036、0.018，每加密一次大致减半；在基础网格上它与 ν=0.5 相当，到最细一层只剩 ν 的 4%。所以观测阶是从下面一步步靠近 1 的：L2 范数 0.59、0.71、0.82、0.90、0.95，最大范数 0.57、0.62、0.80、0.89、0.94（格点页面的 L2 阶是 0.53、0.71、0.82、0.90、0.94）。
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.487, 0.270, 0.140, 0.071, 0.036 and 0.018 level by level, roughly halving with every refinement; on the base mesh it is about as large as ν=0.5, and on the finest it is 4% of ν. That is why the observed order approaches 1 from below, step by step: 0.59, 0.71, 0.82, 0.90 and 0.95 in L2, 0.57, 0.62, 0.80, 0.89 and 0.94 in the max norm (the node-centred page's L2 orders are 0.53, 0.71, 0.82, 0.90 and 0.94).
上面量的是解 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=Σf∈Fwhw，壁面值取 ue，在构造解的离散问题上求值：
(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 sum of the wall-face fluxes J=Σf∈Fwhw, with the wall value ue, evaluated on the discrete manufactured problem:
| 层 ℓ | level ℓ | 单元数 | cells | 只有扩散（β = 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 | 2.94e−3 | — | −5.23e−2 | — | ||||||
| 1 | 96 | 8.18e−4 | 1.85 | −3.23e−2 | 0.70 | ||||||
| 2 | 384 | 2.09e−4 | 1.97 | −1.72e−2 | 0.91 | ||||||
| 3 | 1536 | 5.25e−5 | 1.99 | −8.67e−3 | 0.99 | ||||||
| 4 | 6144 | Jacobi 迭代发散Jacobi diverges | −4.31e−3 | 1.01 | |||||||
| 5 | 24576 | Jacobi 迭代发散Jacobi diverges | −2.14e−3 | 1.01 | |||||||
J 的误差按格式的阶减小：只有扩散时观测阶是 1.85、1.97、1.99，完整方程是 0.70、0.91、0.99、1.01、1.01（格点页面：1.90、1.98、1.99、2.00、2.00；1.05、1.01、1.00、1.00、1.00）。完整方程里，u 的 L2 阶到最细一层是 0.95，还在从下面靠近 1；J 的阶在最细两层是 1.01、1.01，已经是 1。所以对这个目标函数，精度阶就是格式的阶：只有扩散时二阶，完整方程一阶。
The error of J falls at the scheme's order: the observed orders are 1.85, 1.97 and 1.99 with diffusion alone and 0.70, 0.91, 0.99, 1.01 and 1.01 for the full equation (on the node-centred page: 1.90, 1.98, 1.99, 2.00 and 2.00; 1.05, 1.01, 1.00, 1.00 and 1.00). For the full equation the L2 order of u is 0.95 on the finest level, still approaching 1 from below, while the order of J on the two finest levels is 1.01 and 1.01, 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 节在同一族网格上看到两种格式的壁面通量互相靠拢：格心 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 Section 7 shows the two schemes' wall fluxes closing in on each other on this same family: 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. On the node-centred page the multipliers of the strongly imposed rows are not wall field samples; this page imposes the wall weakly and has no such rows. 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.
| 本页（格心） | this page (cell-centred) | 格点页面 | node-centred page | ||
|---|---|---|---|---|---|
| 壁面条件 | wall condition | 弱施加：uw 进入壁面面的通量与模板 | weak: uw enters the wall faces' flux and the stencils | 强施加：壁面节点的行换成 ui−uw | strong: the wall nodes' rows become ui−uw |
| 梯度模板 | gradient stencil | 面邻居，外加壁面面的中点 | face neighbours, plus the wall faces' midpoints | 边邻居 | edge neighbours |
| 目标函数 | objective | 壁面面通量之和 | the sum of the wall-face fluxes | 反力：−Σi∈WRi(u) | reaction: −Σi∈WRi(u) |
| 伴随里的壁面 | the wall in the adjoint | 没有特殊处理 | nothing special | 循环前清零、循环后加回 | zero before the loop, add back after |
| Jacobi 的谱半径 ρ(I − diag(A)−1A)，格点用 Abc | Jacobi spectral radius ρ(I − diag(A)−1A), with Abc on the node-centred page | 0.915 | 0.787 | ||
| 所算最细一层的 L2 阶（只有扩散 / 完整方程） | L2 order on the finest level computed (diffusion alone / full equation) | 2.01 (ℓ = 3) / 0.95 | 2.03 / 0.94 | ||
| 所算最细一层目标函数误差的阶（只有扩散 / 完整方程） | order of the objective's error on the finest level computed (diffusion alone / full equation) | 1.99 (ℓ = 3) / 1.01 | 2.00 / 1.00 |
本项目以离散伴随为主线：从可逐句检查的小模型，逐步走向工业级非结构网格 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.
在同一个多边形区域上用 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.068e-01 | 4.151e-04 | 3.238e+00 |
| 1 | 96 | 4.696e-02 | 4.327e-04 | 3.399e+00 |
| 2 | 384 | 3.277e-02 | 4.390e-04 | 3.392e+00 |
| 3 | 1536 | 2.089e-02 | 4.420e-04 | 3.357e+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∞ 这一列该怎么读。表中第 0 层的 3.238 离参考值 3.292 最近，但这不说明它最准：24 个单元上的误差碰巧互相抵消。从第 1 层起，数值从上方逐层靠近参考值。表只到第 3 层；用本页自己的 Jacobi 迭代（第 10 节的切线求解，每层先到判据、再跑同样多轮）把同一族网格算到第 5 层，得到 3.238、3.399、3.392、3.357、3.329、3.312，前 4 层与表中一致，之后继续逐层下降，第 4 层的 3.329 已经比第 0 层更接近参考值。表中伴随场的 L2 差逐层只缩小 2.27、1.43、1.57 倍，还看不出稳定的收敛率：这几层还在渐近区之前，加密给出的是积分范数下趋向一致的证据，而不是一个量出来的阶。
How to read the dJ/du∞ column. Level 0's 3.238 is the closest to the reference 3.292, but that does not make it the most accurate: on 24 cells the errors happen to cancel. From level 1 on, the values approach the reference from above, level by level. 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.238, 3.399, 3.392, 3.357, 3.329 and 3.312, the first 4 as in the table, and falling further level by level after that; level 4's 3.329 is already closer to the reference than level 0's. The adjoint's L2 difference in the table shrinks only by 2.27, 1.43 and 1.57 per level, which shows no settled rate yet: these levels are pre-asymptotic, and refinement gives evidence of consistency in this integral norm, not a measured order.
结果对不上时应该先查什么（存在性、唯一性、稳定性，边界条件与目标是否匹配，乘子的含义），在附录 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.