方程:一维定常 Burgers 方程
The equation: steady one-dimensional Burgers
从守恒律到定常方程
From a conservation law to the steady equation
一维 Burgers 方程是最简单的非线性守恒律。把 看作沿 方向输运的一个量:一段区间里 的总量,只因两端的通量 和区间内的源项 而改变。写成微分形式,就是
The one-dimensional Burgers equation is the simplest nonlinear conservation law. Think of as a quantity carried along : the amount of in an interval changes only through the flux at its two ends and the source inside it. In differential form,
按乘积法则,: 以自己的值为速度被输运。 正是流体动量方程里的非线性对流项。本页只研究定常解。令 ,得到要在 上离散的方程和它的入口条件:
By the product rule, : is carried along at a speed equal to its own value. The term is exactly the nonlinear convection term of the momentum equation of fluid flow. This page studies steady solutions only. Setting gives the equation to be discretised on and its inlet condition:
我们研究光滑正解分支 上的解。本例不含激波,所以第 3 节的二阶格式使用无斜率限制器的线性重构;间断和激波留到附录 B。
We study solutions on the smooth positive branch . This example has no shocks, so the second-order scheme of Section 3 uses unlimited linear reconstruction; discontinuities and shocks are left to Appendix B.
图 1把这个问题画在一张图上:当前设计,离散解取本页默认的 、一阶(节点值)。下半部分是源项:8 根柱是当前设计的常值 ,虚线是目标源项 。上半部分是解:实线是精确解 ,圆点是本页算出的离散解(第 2 节起讲怎样算),虚线是目标场 。两部分之间的箭头是各处的波速 :全部指向右,并随 增大而变长。入口 给定 ;出口 不给条件。
Figure 1 draws the whole problem: the current design, with the discrete solution at this page's defaults, and first order (nodal values). The lower panel is the source: the 8 bars are the current constant values and the dashed curve is the target source . The upper panel is the solution: the solid curve is the exact solution , the dots are the discrete solution this page computes (Section 2 onwards shows how) and the dashed curve is the target field . The arrows between the panels give the wave speed along the interval: all point right and lengthen as grows. The inlet has the prescribed ; the outlet takes no condition.
波速:信息从入口流向出口
Wave speed: information travels from the inlet to the outlet
通量对状态的导数 就是波速(特征速度): 的一个小扰动以速度 沿 传播。当前解处处 ,所以信息从左向右传播(图 1中间的箭头)。于是左端 是入口,要给一个条件 ;右端 是出口,不给条件,出口值由方程和上游的数据决定。第 4 节说明离散格式怎样处理这两端。
The derivative of the flux with respect to the state, , is the wave speed (the characteristic speed): a small disturbance of travels along with speed . The current solution has everywhere, so information travels from left to right (the arrows in the middle of Figure 1). The left end is therefore the inlet and takes one condition, ; the right end is the outlet and takes none, since the equation and the upstream data determine the outlet value. Section 4 shows how the discretisation treats the two ends.
源项、设计变量与精确解
The source, the design variables and the exact solution
源项在 的 个等长区间上分段取常值 (图 1下半部分的柱);默认 ,这些区间也是第 2 节网格的原始区间。 个源项值与入口值一起构成设计输入 ,改变它们就改变流场:
The source is piecewise constant, with the value on each of equal intervals of (the bars in the lower part of Figure 1); by default , and these intervals are also the primal intervals of the mesh of Section 2. The source values and the inlet value form the design input ; changing them changes the flow:
默认 时, 共有 9 个分量。当前设计取全部 、,方程成为 。它可以直接积分:从 积分到 ,通量的改变量等于源项的累积,再代入常数源项:
With the default , has 9 components. The current design takes for every and , so the equation becomes . It can be integrated directly: from to the change of the flux equals the accumulated source, and with the constant source
取正平方根,因为本页只研究正解分支。出口值 不是给定的,而是由方程和入口数据决定(图 1的实线)。这个闭式解是第 15 节测量网格收敛的参考解,也是第 16 节和附录 A 中连续伴随的出发点。
The positive root is taken because this page studies the positive branch. The outlet value is not prescribed; the equation and the inlet data determine it (the solid curve in Figure 1). This closed form is the reference for the mesh-convergence study of Section 15 and the starting point of the continuous adjoint in Section 16 and Appendix A.
设计的目的在第 7 节给出:调整源项,让解接近目标场 (图 1上半部分的虚线)。目标场由目标源项 和同一个入口值产生;目标函数测量解与目标场的差别,这是一个反设计问题。
Section 7 gives the purpose of the design: adjust the source so that the solution approaches the target field (the dashed curve in the upper part of Figure 1). The target is produced by the target source with the same inlet value; the objective measures the mismatch between the solution and the target, which makes this an inverse-design problem.
这个方程的作用
What this equation is for
在流体力学里,Burgers 方程保留了动量方程中的非线性对流项 ,去掉了压力;它是含有这一项的最简单的方程。没有黏性时,速度大的部分追上前面速度小的部分,波形变陡,最后形成激波(附录 B 讨论激波和它的伴随)。加上黏性项,,变陡与扩散互相平衡,形成有限厚度的光滑过渡层(附录 D.3、D.4 用的就是带黏性的形式)。正因为它简单而又保留了非线性,数值通量、迎风格式和斜率限制器常常先在这个方程上检验。
In fluid dynamics the Burgers equation keeps the nonlinear convection term of the momentum equation and drops the pressure; it is the simplest equation that carries this term. Without viscosity, faster parts of a wave catch up with slower parts ahead of them, the profile steepens and a shock forms (Appendix B treats shocks and their adjoints). With a viscous term, , steepening and diffusion balance in a smooth transition layer of finite width (Appendices D.3 and D.4 use viscous forms of the equation). Because it is simple yet keeps the nonlinearity, numerical fluxes, upwind schemes and slope limiters are commonly tested on it first.
在本课程里,它作为第一个算例,有五个理由:
In this course it is the first example, for five reasons:
- 小。默认网格只有 9 个未知量,每个数组和矩阵都能完整地写在页面上,逐句播放也看得清。
- Small. The default mesh has only 9 unknowns, so every array and matrix fits on the page and the statement players stay readable.
- 非线性。通量的斜率 随状态改变,所以线性化(Jacobian)取决于在哪个状态上求导:前向与伴随都在收敛的原始解上计算。
- Nonlinear. The slope of the flux, , changes with the state, so the linearisation (the Jacobian) depends on the state it is taken at: the tangent and the adjoint are both evaluated at the converged primal solution.
- 有精确解可比。闭式解 用来测量离散误差(第 15 节);连续伴随也能写出来(第 16 节、附录 A),用来和离散伴随对照。
- Exact answers to compare with. The closed form measures the discretisation error (Section 15), and the continuous adjoint can be written down too (Section 16, Appendix A) for comparison with the discrete one.
- 设计变量在源项里。改变设计不会移动网格,所以不需要对网格坐标求导;几何导数留给二维页面。
- Design variables in the source. Changing the design does not move the mesh, so no derivatives with respect to mesh coordinates are needed; geometric derivatives are left to the 2-D pages.
- 同一个面循环。本页的面循环(收集、重构、算通量、分发)在 Euler、N–S 和 RANS 页面上反复出现,只是每个面上的计算更多;黏性页面还要先算一遍梯度。
- The same face loop. This page’s face loop (gather, reconstruct, evaluate the flux, scatter) reappears on the Euler, N–S and RANS pages, with more work on each face; the viscous pages first compute gradients.
图 2画出通量 。当前解的状态从入口的 增加到出口的 ,占据曲线上加粗的一段。两端切线的斜率就是两处的波速:入口处为 ,出口处为 。斜率随状态改变,这就是非线性:在一个状态上求出的导数,不能拿到另一个状态上用。所以第 3 节数值通量的偏导 和第 6 节的 Jacobian 都要在当前状态上计算,前向与伴随则在收敛的原始解上计算。竖直方向,两端通量之差 恰好等于源项的积分 :定常时,流出的通量减去流入的通量,等于区间里的源项。第 3 节就从这个守恒关系出发,把它写到每一个控制体上。
Figure 2 plots the flux . The state of the current solution rises from at the inlet to at the outlet and covers the thick part of the curve. The slopes of the tangents at the two ends are the wave speeds there: at the inlet and at the outlet. The slope changes with the state, and that is the nonlinearity: a derivative computed at one state cannot be reused at another. The flux derivatives of Section 3 and the Jacobian of Section 6 are therefore evaluated at the current state, and the tangent and the adjoint at the converged primal solution. Vertically, the difference of the end fluxes, , equals the integrated source : in a steady state the flux out minus the flux in equals the source inside. Section 3 starts from this conservation statement and writes it for every control volume.
图 3把本页放在整个网站的路线上。每往右一步主要增加一类新内容:拟一维 Euler 把一个方程换成质量、动量、能量三个方程,并引入压力和变截面积;二维 Euler 增加第二个空间方向,面有了法向,网格坐标也成为设计变量;二维层流 N–S 加入黏性应力和热传导,面上还要用到梯度;二维 RANS–SA 再加一个湍流模型方程。每点的未知量从 1 个增加到 5 个,但每一页的残差仍由面循环组装,前向与伴随也沿同一个循环求导。另有两页与本页相邻:一维层流 N–S 页留在一维,在 Euler 方程上加入黏性和热传导;二维标量对流页是本页在二维的对应页面。
Figure 3 places this page on the route through the site. Each step to the right mainly adds one kind of ingredient: quasi-1-D Euler replaces the single equation by three (mass, momentum and energy) and brings in the pressure and a varying cross-section; 2-D Euler adds a second space direction, so faces have normals and the mesh coordinates become design variables; 2-D laminar N–S adds viscous stresses and heat conduction, which need gradients at the faces; 2-D RANS–SA adds one turbulence-model equation. The unknowns per point grow from 1 to 5, yet every page still assembles its residual with face loops, and the tangent and the adjoint differentiate the same loops. Two more pages sit next to this one: the 1-D laminar N–S page stays in one dimension and adds viscosity and heat conduction to the Euler equations, and the 2-D scalar advection page is this page's counterpart in two dimensions.
本页要回答的问题
The question this page answers
假设你已经学过微积分、矩阵乘法、Newton 法和线性方程的迭代求解。本页要连接这些知识:先把一个微分方程写成有限个非线性代数方程,再求“输入改变一点,解和目标函数会改变多少”。你不需要先学自动微分或最优化;前向与伴随都从链式法则和矩阵转置开始。
This page assumes calculus, matrix multiplication, Newton’s method and iterative linear solvers. We connect them by turning a differential equation into finitely many nonlinear equations, then asking how a small input change affects the solution and an objective. Automatic differentiation and optimisation are not prerequisites: forward and adjoint calculations start with the chain rule and matrix transposition.
本页只回答一个问题:设计输入 改变一点,目标 改变多少;页首的全页路线图把后面各节串成一条线。按页首默认设置(一阶、GS 固定点、每次一轮扫描)实际计数,9 次前向求解共用 210 次 GS 扫描,1 次伴随求解只用 24 次。
This page answers one question: how much does the objective change when the design input changes a little? The page roadmap at the top strings the later sections together. Counted with the default settings at the top of the page (first order, GS fixed point, one sweep per application), the 9 forward solves take 210 GS sweeps in total and the single adjoint solve only 24.
要在计算机上回答这个问题,先要把连续方程换成有限个未知量的代数方程。第 2 节先建立网格、决定未知量放在哪里;第 3 节再把守恒律写成面上的通量差。
To answer it on a computer, the continuous equation must first become algebraic equations in finitely many unknowns. Section 2 builds the mesh and decides where the unknowns live; Section 3 then writes the conservation law as flux differences across faces.
网格:未知量放在哪里
The mesh: where the unknowns live
第 1 节的方程对区间上的每一个 都成立,计算机却只能保存有限个数。本节选网格:决定未知量放在哪里、数组怎样编号,再写下本页用到的符号和数据。
The equation of Section 1 holds at every of the interval, but a computer stores finitely many numbers. This section chooses the mesh: it decides where the unknowns live and how the arrays are numbered, and fixes the notation and data used on this page.
先把物理问题、网格和数组分开
Separate the physical problem, mesh and arrays
用 表示原始区间数, 表示网格宽度, 表示未知量个数。数组从 0 开始编号。格心格式有 ;格点格式有 。网格位置 是已知几何数据,未知量 是要解出的数值,二者不能混淆。
Let be the number of primal intervals, the spacing and the number of unknowns. Arrays use zero-based indices. Cells have unknowns; nodes have . Coordinates are known geometric data; values are numbers to solve for.
格点格式先在相邻节点的中点切开区间。围绕一个内部节点的两半段合成它的对偶控制体;端点只有半个控制体。源项仍存储在原始区间上,所以一个内部节点的源项积分需要左右两个区间的信息。入口节点的值已给定,我们用一条等式约束代替它的平衡方程。
For a node scheme, split each primal interval at its midpoint. The adjacent half-intervals form the dual volume around an interior node; endpoint volumes have half the width. Sources remain stored on primal intervals, so an interior node receives source integrals from two intervals. The prescribed inlet node uses an equality constraint in place of its balance equation.
本页的符号和数据
Notation and data used on this page
是状态向量, 是设计输入, 是残差, 是一个标量目标。按伴随优化文献的通常写法,设计变量记作 。求导采用自动微分(AD)的记号:程序中任一变量 的前向(tangent)导数记作 ,反向(adjoint)累加到 上的权重记作 。设计方向是种子 ; 是状态的切向量,也就是前向求解得到的数组; 是检验转置乘积时任意给定的残差种子 , 是由伴随方程求出的特定种子。 和 的用途不同。
The state is , the design input is , the residual is , and the objective is a scalar. As is usual in the adjoint-optimization literature, the design variables are written . Derivatives use the notation of automatic differentiation (AD): the forward (tangent) derivative of any program variable is , and the weight accumulated into in reverse (adjoint) mode is . A design direction is the seed ; is the state tangent, the array returned by a tangent solve; is an arbitrary residual seed for testing a transpose product, and is the particular seed obtained by solving the adjoint equation. These two seeds serve different purposes.
当前计算网格
The current computational mesh
未知量是 处的节点值,。残差在对偶控制体上积分;端点控制体宽度为 ,内部为 。入口第 0 行稍后替换为约束。
Unknowns are point values at , . Residuals are integrated over dual volumes: width at the endpoints and inside. Row 0 will be replaced by the inlet constraint.
页首的一阶/二阶按钮控制整页计算。求解器选择独立于空间精度。所有实时示例使用当前网格、边界处理与设计变量;精度研究使用固定的网格序列,并同时画出两个阶数。
The first-/second-order buttons control calculations throughout this page. Solver choice is independent of spatial order. Live examples share the selected mesh, boundary treatment and design; the accuracy study uses a fixed mesh sequence and displays both orders.
网格与面:计算前准备的数据
Mesh and faces: data prepared before calculation
h = 1 / N; M = N + 1
x[i] = i*h # i = 0,...,N
W = [h/2, h, ..., h, h/2]
faces = [(i, i+1, (i+1/2)*h) for i = 0,...,N-1]
faces += [(N, -1, 1)] # physical outflow
# No exterior inlet face: row 0 is replaced by u[0] - u_in.伪代码约定:一个面记为 ; 表示域外,不是合法数组下标。只有存在对应控制体时才向它累加。 是控制体宽度;后面的目标积分也使用它。
Pseudocode convention: a face stores . The index means exterior, not a valid array index. Accumulate only when the corresponding volume exists. The width also supplies the objective quadrature weight.
网格定下了控制体和它们的面。第 3 节把第 1 节的守恒律写到控制体上:控制体里源项的总量,等于从它两端的面流出的净通量;所以先要会算一个面上的通量。
The mesh fixes the control volumes and their faces. Section 3 writes the conservation law of Section 1 on the control volumes: the total source in a volume equals the net flux out through the faces at its two ends, so the first task is the flux on one face.
面上的通量:守恒与重构
The flux on a face: conservation and reconstruction
本节实验:改变一个未知量,哪些面值会变化?先预测重构邻域,再用图验证。
Experiment: Which traces change when one unknown changes? Predict the stencil, then inspect the diagram.
第 2 节把区间分成了控制体。本节把守恒律写成控制体上的通量差,用数值通量把面两侧的状态变成一个面通量,再说明二阶格式怎样把状态重构到面上。
Section 2 divided the interval into control volumes. This section writes the conservation law as a flux difference over a control volume, turns the states on the two sides of a face into one face flux with a numerical flux, and then shows how the second-order scheme reconstructs states at the faces.
从守恒到通量差
From conservation to a flux difference
物理通量就是第 1 节的 。对控制体 积分,用微积分基本定理,导数积分变为两个端点的通量差。这一步是有限体积法的起点:控制体内的源项总量,必须等于净流出通量。
The physical flux is the of Section 1. Integrate over a volume and use the fundamental theorem of calculus: the derivative becomes the difference of endpoint fluxes. This is the starting point of finite volumes: the integrated source equals the net outward flux.
离散后,相邻控制体可能给同一个面提供不同的近似值 ,所以需要一个数值通量 把这两个值变成一个面通量。这里的 是状态,不是上一式的积分端点。通量应满足一致性 。共享同一个 ,还让内部面的贡献在全域求和时严格相消。
After discretisation, neighbouring volumes can supply different states at the same face. A numerical flux turns these into one shared flux. Here denote states rather than integration endpoints. Consistency requires . Sharing the same flux also makes interior contributions cancel in a global sum.
本例固定耗散系数为 。项 抑制相邻状态的跳跃;它也会影响离散误差。第 1 节的精确解 就是后面测量这种误差的参考解。
The dissipation coefficient is fixed at . The term damps state jumps and affects discretisation error. The exact solution of Section 1 is the reference against which that error is measured later.
从一阶常值到二阶重构
From constant states to linear reconstruction
一阶格式在控制体内取常值。二阶内部斜率用中心差分,端点斜率用单边差分。虽然端点斜率本身只有一阶精度,但乘以 的外推距离后,面值仍有 精度;实际全域阶数由网格加密验证。
First order uses piecewise constant values. Second order uses centred interior slopes and one-sided endpoint slopes. The endpoint slope is only first-order accurate, but multiplication by an extrapolation distance gives an face value; mesh refinement checks the resulting global order.
其中 是未知量个数。通量的两个偏导为 、。重构之后仍使用同一个通量;二阶精度来自面值重构。下图可选择一个面,查看两侧斜率与重构值。
Here is the number of unknowns. Flux derivatives are and . Both orders use the same numerical flux: reconstruction changes the face values. Select a face below to inspect the slopes and reconstructed states.
把重构写成可以求导的加权和
Write reconstruction as differentiable weighted sums
二阶不是更换 Burgers 通量,而是用相邻数据估计面上的状态。以内部左迹为例,把中心斜率代入 ,就得到三个非零权重。右迹同理;合起来一个面依赖四个未知量。这也解释了为什么二阶 Jacobian 比一阶更宽。
Second order estimates face states using neighbouring data while retaining the Burgers flux. Substitute the centred slope into to obtain three nonzero weights. The right trace works similarly; together a face depends on four unknowns. This explains the wider second-order Jacobian.
重构伪代码:一阶与二阶共用入口
Reconstruction pseudocode: one interface for both orders
trace(i, xf, u, order):
if order == 1: return u[i]
if i == 0: slope = (u[1] - u[0]) / h
else if i == M-1: slope = (u[M-1] - u[M-2]) / h
else: slope = (u[i+1] - u[i-1]) / (2*h)
return u[i] + (xf - x[i]) * slope读图练习:先选一阶,再选二阶,保持面编号不变。观察左右迹、模板点与通量怎样变化。线性函数应被二阶重构精确再现;这说明局部重构正确,但全局二阶还需要边界、源项积分和网格加密结果共同支持。
Try keeping the same face selected while switching order. Observe the traces, stencil points and flux. Second-order reconstruction reproduces a linear function exactly. This checks local reconstruction; global second order also depends on boundaries, source integration and mesh refinement.
一个面上的导数:通量与重构
Derivatives on one face: the flux and the reconstruction
第 8、9 节的前向与伴随只用到两类局部导数:数值通量对面两侧状态的偏导,以及面值(迹)对未知量的偏导。这里把它们一次推导清楚。数值通量是两侧状态的多项式,逐项求导即可:
The tangent and the adjoint of Sections 8 and 9 use only two kinds of local derivative: those of the numerical flux with respect to the states on the two sides of a face, and those of the face values (the traces) with respect to the unknowns. Derive them once here. The numerical flux is a polynomial in its two states, so it can be differentiated term by term:
出口面(以及“边界对照”里的迎风面)用物理通量 ,它只依赖左侧状态:,。有两条性质可以检查这些公式。一致性 对每个 都成立,沿 求导得 :两个偏导把波速分成上游与下游两份。在本页的状态范围 上, 在 与 之间, 在 与 之间:通量随上游状态增大、随下游状态减小。只要 ,这种单调性就成立。
The outlet face (and the upwind face of the boundary comparison) uses the physical flux , which depends on the left state only: and . Two properties check these formulas. Consistency holds for every ; differentiating it along gives : the two partial derivatives split the wave speed into an upstream and a downstream share. Over this page's states , lies between and and between and : the flux increases with the upstream state and decreases with the downstream one. This monotonicity holds as long as .
面值的导数更简单。不用斜率限制器时,每个迹都是未知量的固定线性组合 ,所以它对 的偏导就是权重本身,与 无关:一阶时所在控制体的权重为 ;二阶内部面是上式的 ,第一个和最后一个内部面用单边斜率,权重为 ;格点出口节点就在 ,它的迹权重为 。权重不随状态变化,所以前向沿用同一组权重,伴随用它们的转置,并用 += 累加:
The traces are simpler still. Without a limiter each trace is a fixed linear combination of the unknowns, , so its partial derivative with respect to is the weight itself and does not depend on : weight on its own volume at first order; at second order for an interior trace as above, at the first and the last interior faces, where the slope is one-sided, and the outlet node lies at , so its trace has the weight . Since the weights do not change with the state, the tangent reuses the same weights and the adjoint their transpose, accumulated with +=:
本页不用斜率限制器。若采用受限重构 (附录 D.2),限制器 随状态变化,切向量多出一项:。权重于是依赖 ,而 minmod 在切换分支处不可微(附录 D.1)。
This page uses no limiter. A limited reconstruction (Appendix D.2) has a limiter that changes with the state, so its tangent gains a term, : the weights then depend on , and where minmod switches branches it is not differentiable (Appendix D.1).
图 6把这两类导数画出来:左边是 、 随状态的变化;右边是本页默认网格上每个面的通量对每个未知量的导数 。一阶时每个内部面只连两个未知量,二阶时连四个(第一个和最后一个内部面只连三个);第 6 节把这些行按“左加右减”组合成 Jacobian。
Figure 6 draws both kinds: on the left and against the state; on the right, for this page's default mesh, the derivative of every face flux with respect to every unknown, . At first order an interior face touches two unknowns, at second order four (three at the first and the last interior face); Section 6 combines these rows, plus on the left and minus on the right, into the Jacobian.
到这里,每个内部对偶面的通量都能由两侧的重构状态算出。区间的两端不同:入口节点 的值已经给定,它的行要换成约束;出口节点 的对偶控制体右端就是区间端点,那里只用物理通量。第 4 节写出这两端的处理。
Every interior dual face now gets its flux from the reconstructed states on its two sides. The two ends of the interval differ: the value at the inlet node is prescribed, so its row becomes a constraint, and the dual volume of the outlet node ends at the end of the interval, where only the physical flux is used. Section 4 writes how the two ends are treated.
边界条件:入口与出口
Boundary conditions: inlet and outlet
为什么只给一个入口条件
Why prescribe only one inlet condition?
第 1 节从波速的方向说明了只给入口条件的原因。换一个角度看,定常方程在正解分支上是一阶常微分方程。给定 后,积分就确定沿程的 ;出口值应由内部解决定。再任意指定一个出口值,通常会过约束。离散时仍须提供出口通量,但“计算出口通量”与“另加出口值约束”是两件事。
Section 1 explained the single inlet condition by the direction of the wave speed. Seen another way, the steady equation on the positive branch is first order. Once is specified, integration determines downstream. An arbitrary additional outlet value would generally overconstrain it. The discretisation still needs an outlet flux, but evaluating that flux is different from imposing an extra outlet-value constraint.
入口保留强约束 。第一内部对偶面位于 ,二阶时必须重构到这个位置;不能仍用 充当面值。出口节点已在 ,直接使用物理通量 。
The inlet retains the strong constraint . The first internal dual face is at ; second order must reconstruct to that location rather than use as the face value. The outlet node is already at , so its physical flux is .
图 8用本页默认网格(,一阶)的编号和收敛值画出第 0 行的处理顺序。面循环照常把第一个内部面 的通量加到第 0 行;循环结束后,第 0 行整行替换为 ,循环写入的值被丢弃。伴随按相反顺序做转置:先令被替换行的残差种子 (原来的 另外留着),再做转置面循环,最后补上约束行的反向语句 。严格按原始语句的反序,这条语句是第一条;+= 可以交换次序,只要 在面循环读取它之前已经清零,结果相同,图中按上面的顺序画。最后一步必须累加,不能覆盖面循环已经传回 的量。
Figure 8 draws the order of operations for row 0 with the indices and converged values of the default mesh (, first order). The face loop adds the flux of the first interior face to row 0 as usual; after the loop, row 0 is replaced by and the value the loop wrote is discarded. The adjoint applies the transposed steps in reverse order: first set the residual seed of the replaced row to (the original is kept aside), run the transposed face loop, and finally add the reverse statement of the constraint row, . In strict reverse order of the primal statements this one comes first; += commutes, so the result is the same as long as is cleared before the face loop reads it, and the figure draws the order given here. The last step must accumulate; it must not overwrite what the face loop has already returned to .
“边界对照”允许改用迎风通量 :格心改入口面,格点改第一内部对偶面。它是另一种明确标记的离散边界闭合,用于附录 A 的比较。格点二阶迎风面依赖重构斜率,因此它仍依赖相邻未知量。
The boundary comparison replaces one flux by : the inlet face for cells, or the first internal dual face for nodes. This is a distinct, explicitly labelled boundary closure, compared with the default one in Appendix A. At second order, the node-based upwind face depends on its slope and therefore still on neighbouring unknowns.
对状态求导时,约束行 中的 固定,只剩 ;对入口设计变量求导时,。这两个导数分别进入 和 ,不能混在一个状态模板里。边界对照改变的是离散方程本身,因此必须同时改变原始、前向、伴随和目标场。
When differentiating with respect to the state, in the constraint row is fixed, leaving ; with respect to the inlet design variable, . These paths enter and separately. Changing the boundary comparison changes the discretisation and must be reflected in the primal, forward, adjoint and target field.
边界行怎样进入前向与伴随
How the boundary rows enter the tangent and the adjoint
格点页的两条边界行形式不同:第 0 行被约束整行替换;最后一行属于出口节点,它只有半个控制体,左侧一个内部面带负号,出口面(物理通量)带正号:
The two boundary rows of the node page differ in kind: row 0 is replaced as a whole by the constraint; the last row belongs to the outlet node, whose volume is a half, with one interior face on its left (minus) and the outlet face (plus, the physical flux):
前向:约束行对 与 的偏导是 与 ;出口节点在 ,它的迹就是 本身:
Forward: the constraint row has the partial derivatives and with respect to and ; the outlet node lies at , so its trace is itself:
伴随:约束行在原始代码里最后写入,所以严格按反序它的反向语句最先执行,并在用完 后把它清零,此前写入第 0 行的通量与源项不再得到权重(第 9 节的播放器与图 8 把 放在面循环之后:+= 可以交换次序,只要 在面循环读取它之前已经清零,结果相同);出口行照常累加,它含出口面和最后一个内部面,这个内部面的两个迹是 与 :
Adjoint: the constraint row is written last in the primal code, so in strict reverse order its statement runs first and clears once it has used it; the flux and the source written into row 0 before receive no weight (the players of Section 9 and Figure 8 place after the face loop instead: += commutes, so the result is the same as long as is cleared before the face loop reads it). The outlet row accumulates as usual; it holds the outlet face and the last interior face, and the two traces of that interior face are and :
在默认网格()的收敛解上,这两行的非零元素如下表,它们就是第 6 节 Jacobian 的第 0 行、最后一行与第 10 节 的对应元素。上面的语句可以直接核对:取 只得到 ,一阶时取 得到 、,正是表中的 、 与 。第 0 行对源项的偏导为零:约束行不含源项。二阶时出口行多出 。
At the converged solution on the default mesh () the nonzero entries of the two rows are listed below: they are row 0 and the last row of the Jacobian of Section 6 and the matching entries of of Section 10. The statements above can be checked directly: gives only , and at first order gives and , which are the entries , and of the table. Row 0 has no source derivative: the constraint row carries no source. At second order the outlet row gains .
| 元素Entry | 一阶First order | 二阶Second order |
|---|---|---|
| — | ||
残差:收集、重构、计算、分发
Residual: gather, reconstruct, compute, scatter
先固定离散问题与符号
Fix the discrete problem and its notation
本页每个节点存储 1 个未知量。以下用 表示把所有位置按代码顺序排成的列向量(标量问题每个位置只有一个分量)。它是一个有限维向量,不是连续函数。物理方程中的空间导数只用于说明模型;求离散伴随时,求导对象是本页实际装配的有限维函数。
Each node stores 1 unknown. Below, stacks all locations in code order (a scalar problem has one component per location). It is a finite-dimensional vector, not a continuous function. Spatial derivatives in the physical equations explain the model; discrete differentiation acts on the actual finite-dimensional functions assembled on this page.
是本页全部设计变量排成的向量:每个原始区间上的常源项 与入口值 。第 6 节按空间阶数把精确 Jacobian 记作 ,这里的 就是当前阶数的 ;第 7 节的 就是这里的状态梯度;第 10 节的 是 作用在设计种子 上的结果。
stacks all design variables of this page: the constant source on each original interval and the inlet value . Section 6 writes the exact Jacobian of each spatial order as ; the here is the of the current order. The of Section 7 is the state gradient here, and the of Section 10 is applied to the design seed .
这里 、 都是列向量,它们的转置才是偏导数行。求导采用自动微分(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.
残差是“还没有平衡的量”
A residual measures the remaining imbalance
把所有离散平衡式的右侧移到左侧,得到向量 。一个分量对应一个方程; 表示该控制体平衡。求解器的任务是找使全部分量同时为零的 。计算残差本身只是在给定 后做一次函数求值,不会自动求出解。
Move each discrete balance to the left to form . One component represents one equation; means balance in that volume. The solver seeks a state making all components zero simultaneously. Evaluating the residual at a supplied state is only a function evaluation, not a solution of the equations.
面循环比逐行重复计算通量更直接:先把残差数组清零,再遍历每个面,重构状态、计算通量、向两侧累加,最后扣除源项。内部通量在两个方程中重复出现,却只计算一次。这里存储的是积分残差,没有除以控制体宽度 ;后面的 Jacobian、目标权重和伴随符号都按这一约定。
A face loop evaluates each shared flux once: zero the residual, visit each face, reconstruct states, compute the flux, add it to its neighbours, then subtract sources. An interior flux appears in two equations but is evaluated once. We store integrated residuals without dividing by volume width ; all Jacobians, objective weights and adjoint conventions follow this choice.
三个矩阵符号分别做什么
What the three matrix symbols do
的一行存储一个面的重构权重;一阶每行通常只有一个 1,二阶每行通常有三个非零数。已知的入口常数另存为 。 的一列对应一个面,在左行放 、右行放 。所以矩阵写法只是把刚才的收集与分发打包,并没有增加新的数值方法。
store one face’s reconstruction weights per row: typically one unit entry at first order and three nonzero entries at second order. Prescribed inlet constants are stored separately in . A column of represents one face, with in the left row and in the right. The matrix notation packages gather and scatter; it introduces no new numerical method.
这些矩阵式首先描述未经入口行替换的通量部分。出口只用左迹的物理通量;格点的第 0 行最后单独施加约束。播放器中的权重列表就是这些稀疏矩阵的非零项。
These matrix expressions first describe the flux contribution before inlet row replacement. The outlet uses the physical flux of its left trace, and the node scheme applies its row-0 constraint separately. The player’s weight lists are the nonzero entries of these sparse maps.
设计变量为每个原始区间上的常源项 和入口值。源项先积分到控制体,面通量按左正右负分发。符号约定固定为“流出减流入减源项”。
Design variables are the constant source values on the primal intervals and the inlet value. Integrate sources over control volumes, then scatter each face flux with a positive sign to its left and a negative sign to its right. The residual convention is outflow minus inflow minus source.
格点源项公式中不存在的区间项取零。将重构写为线性映射 ,将面到控制体的有向分发写为 ,可以把两个空间阶数统一为下面的计算链。边界常数和入口行替换需要额外处理。
Missing interval terms in the node-based source formula are zero. Linear reconstruction maps and the oriented face-to-volume scatter unify both spatial orders. Boundary constants and the inlet row replacement are handled separately.
怎样使用逐句播放器:高亮行是刚执行的语句,变量表只显示已经算出的量,下层网格是执行后的累积数组。“下一步”前进一句;同一面内先得到左右迹,再算通量,最后分别更新左右残差。源项和入口约束也有独立步骤。符号 只表示域外,访问数组之前会被排除。
Reading the statement player: the highlighted line has just executed, the variable table contains values already computed, and the lower mesh shows the accumulated array afterwards. Next advances one statement: left trace, right trace, flux, left update, right update. Sources and the inlet constraint have separate steps. Exterior index is excluded before array access.
Jacobian 与 Gauss–Seidel 迭代
The Jacobian and Gauss–Seidel iteration
从 Newton 法得到线性修正方程
From Newton’s method to a linear correction
你熟悉的一元 Newton 法用切线逼近函数。对向量残差也是同样的思路:在当前状态附近作一阶 Taylor 展开,让线性近似的残差为零。矩阵 表示第 个未知量增加一点,会让第 个残差改变多少。修正量记作 :它是求解器每一步的有限增量,不是求导记号。
Scalar Newton uses a tangent line to approximate a function. For a residual vector, use a first-order Taylor expansion at the current state and set the linear approximation to zero. The entry measures how residual responds to a small change in unknown . The correction is written : a finite solver increment, not a derivative symbol.
精确 Newton 会使用 并充分求解这个线性方程。本页为了展示实际求解器中常用的低阶近似,改用 ,再用有限次 GS 近似解修正方程。每次更新后仍重新计算完整的 。因此必须分清三层:空间格式决定 ;线性化决定 ;近似矩阵和 GS 决定修正怎样得到。
Exact Newton uses and solves that linear system accurately. To illustrate a common low-order approximation, this page uses and a finite number of GS sweeps. Every update is checked with the full residual . Distinguish three levels: spatial discretisation defines , differentiation defines , and the approximate matrix plus GS determines the correction.
GS 中“立即使用新值”是什么意思
What does GS mean by using new values immediately?
第一轮先算 ,随后立即用这个新值算 。第二轮得到 、,逐步接近精确解 。如果第二行仍用旧的 ,就不是这次 GS 扫描。下方网格用绿色标出已经更新的位置。
The first sweep computes , then immediately uses it in . The next sweep gives and , approaching . Using the old in the second row would not be this GS sweep. The mesh below marks updated positions in green.
一个内部残差行怎样组成 Jacobian
How one interior residual row forms the Jacobian
先取一阶内部行。它包含右面的 和左面的 ;分别求导并保持符号,就得到三对角的三个系数。二阶时每个面还依赖远邻,通过重构权重再乘一次,就扩展到五对角。
For a first-order interior row, differentiate the right-face and the left-face with their signs intact. This gives the three tridiagonal coefficients. At second order each face also depends on farther neighbours; multiplication by reconstruction weights widens the row to five diagonals.
这个公式只针对两个相邻面都用内部数值通量的平衡行。入口、出口和约束行必须由各自实际表达式求导,不能直接照抄内部公式。
This formula applies to a balance row whose two faces both use the interior numerical flux. Inlet, outlet and constraint rows must be differentiated from their actual expressions.
精确 Jacobian 是当前空间格式的残差导数 。一阶内部通常是三对角,二阶重构扩大为五对角;边界单边模板也必须求导。求解时只用在当前状态上计算的一阶 Jacobian 构造修正步骤。
The exact Jacobian is for the selected spatial order. Interior first-order rows are generally tridiagonal; reconstruction widens second-order rows to five diagonals. One-sided boundary stencils must also be differentiated. The correction uses the first-order Jacobian evaluated at the current state.
一轮 GS 扫描怎样更新未知量
How one GS sweep updates the unknowns
把 分为对角、严格下三角和严格上三角部分。内层从零修正开始做固定次数 GS 扫描;一轮扫描中,已经经过的行使用新值。默认一次扫描,属于不精确的线性求解;因此需要继续外层残差校正。
Split into diagonal, strictly lower and strictly upper parts. The inner solve starts from zero correction and performs a fixed number of GS sweeps. Already visited rows use new values. The default single sweep is an inexact linear solve, so outer residual corrections continue.
求解用的矩阵与要求解的方程。一阶 Jacobian 给出修正方向,完整残差 决定何时收敛。二阶求解因此仍保留全部重构项。若修正使状态不再为正,就减小步长以留在当前解支。
The correction matrix and the equation being solved. The first-order Jacobian supplies a correction; the full residual determines convergence. A second-order solve therefore retains every reconstruction term. Updates are damped when needed to stay on the positive solution branch.
下面先看两个矩阵的非零结构,再逐行播放一次 GS 扫描,最后观察解和残差怎样随迭代变化。
Below, compare the two matrix stencils, play one GS sweep row by row, then watch the solution and residual change over the iterations.
原始求解:外层更新与内层 GS
Primal solve: outer updates and inner GS
u = positive_initial_state()
repeat:
r = residual(u, D, selected_order)
if norm(r) <= tolerance: return u
P = first_order_jacobian(u)
delta = zeros(M)
repeat fixed_number_of_sweeps:
for i = 0,...,M-1:
delta[i] = (-r[i] - Σ(j≠i) P[i,j]*delta[j]) / P[i,i]
alpha = relaxation
while min(u + alpha*delta) <= 0.05 and alpha > 1e-8: alpha = alpha/2
u = u + alpha*delta
# At the iteration limit, report non-convergence.数组 delta 存放修正量 。GS 行更新中的同一个数组同时保存“前面行的新值”和“后面行的旧值”。本页每次外层修正把 delta 重新置零;同一次修正内的第二轮 GS 则保留第一轮结果。默认一轮 GS 不等于精确求逆。
The array delta holds the correction . Within a GS row update, one array holds new values for visited rows and old values for unvisited rows. Each outer correction starts delta at zero; a second inner sweep retains the first sweep’s result. The default single sweep is not an exact inverse.
下面的动画把原始 GS 缺陷校正画成曲线。先逐行播放第一次修正:正向 GS 扫描从第 0 行走到最后一行,已经更新的一段先离开初值;之后每帧是一次外层迭代。虚线是收敛解 ,右边是残差 的历史,下方的柱是各控制体的误差(对数刻度)。右图竖虚线是页面求解器按停止准则结束的位置;动画继续迭代,残差停在舍入平台上,说明迭代确实收敛到了机器精度。比较入口和出口两侧的柱哪边先变短:原始信息沿特征速度 从入口流向出口。
The animation below draws the primal GS defect correction as a curve. The first correction is played row by row: the forward GS sweep runs from row 0 to the last row, and the updated part leaves the initial value first. After that, each frame is one outer iteration. The dashed curve is the converged solution , the right panel is the history of , and the bars below show the error in each control volume on a log scale. The vertical dashed line on the right marks where the page's solver stops by its criterion; the animation keeps iterating until the residual rests on its round-off floor, which shows that the iteration has converged to machine precision. Compare which side of the bars shrinks first: primal information travels with the characteristic speed from the inlet to the outlet.
符号提示:线性系统中的 (按空间阶数写作 )是残差 Jacobian; 是一阶近似矩阵; 是从零修正开始、固定 轮 GS 扫描构成的线性映射,一般不等于 。
Notation: in a linear system (written for spatial order ) is the residual Jacobian; is the first-order approximate matrix; is the linear map of GS sweeps started from a zero correction, generally not .
目标函数:希望解变成什么样
The objective: the field we want to match
把整条曲线的差异变成一个数
Turn a field mismatch into one number
如果要调整源项,让当前解接近给定目标场,就需要一个衡量差异的标量。平方差不受误差正负抵消影响,控制体宽度 使求和近似空间积分。前面的 只为了求导时消去平方带来的系数 2。
To adjust sources toward a target field, measure the mismatch by a scalar. Squaring prevents positive and negative errors from cancelling; volume weights make the sum approximate an integral. The factor cancels the factor 2 when differentiating the square.
这里 只描述目标函数对状态的直接敏感度,还不是对源项的梯度;目标场 冻结,没有切向量。源项改变后,所有状态会通过平衡方程联动;我们还需要计算 ,或用伴随把这层依赖消掉。
The vector describes direct sensitivity to the state, not yet to the sources; the target field is frozen and has no tangent. A source change affects the state through the balance equations. We must either compute or use an adjoint to eliminate that dependence.
目标是让解接近目标源项 产生的流场。每种离散格式先求出自己的目标场 ,然后在设计求导时保持它不变。切换阶数或边界格式会重建目标场;当前设计改变时它不重新求解。
The inverse objective matches the field generated by . Each discretisation first computes its own target field , which is then frozen during design differentiation. Changing spatial order or boundary closure rebuilds this target; changing the current design does not.
因此,在目标设计处,离散 和梯度同时为零。第 15 节的精度比较另用统一连续参考,防止把各网格自己的目标当作误差基准。
Thus the discrete objective and gradient vanish at the target design. Section 15 uses a shared continuous reference for its accuracy comparisons, rather than treating each grid’s discrete target as its own error reference.
目标函数与状态梯度
Objective and state gradient
target = solve_primal(target_source, prescribed_inlet) # compute once
J = 0
for i = 0,...,M-1:
mismatch = u[i] - target[i]
J += 0.5 * W[i] * mismatch*mismatch
g[i] = W[i] * mismatch
# Reuse this frozen target for every design change.从离散目标得到伴随右端
From a discrete objective to the adjoint right-hand side
不能只对连续积分写一个形式导数就结束。先确定代码实际使用的求积点、权重、重构状态及边界值,再对这个标量表达式求导。一般的求积和及其链式法则是:
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:
平方跟踪目标是它的特例:, 是被跟踪量减去它的目标值,于是 ; 在被跟踪量的切向量之外,只多出目标值本身的切向量(目标冻结时为零)。
A squared tracking objective is the special case , where is the tracked quantity minus its target; then , and adds to the tangent of the tracked quantity only the tangent of the target itself (zero for a frozen target).
本页第 7 节的跟踪目标 正是平方跟踪型:求积点就是存储点,,;目标场 冻结,所以 。它不显式依赖 ,所以 。
The tracking objective of Section 7 is exactly of this form: the quadrature points are the stored points, and ; the target field is frozen, so . It has no explicit dependence on , so .
求 时固定所有设计变量;求 时固定全部离散状态槽位。本页网格不随设计变化,采样点固定,所以没有位置依赖项。目标中的体积权重已经进入 ,不应在伴随右端再乘一次体积。
Hold every design variable fixed for , and every stored state entry fixed for . The mesh of this page does not change with the design, so the samples are fixed and there is no location dependence. 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, .
前向:沿着面循环传播切向量
Forward: propagate tangents through the face loop
本节实验:一次 tangent 计算给出残差对指定方向的变化,不是已经求出了状态对设计的敏感度。比较第 10 节的线性求解。
Experiment: One tangent application gives the residual change along a supplied direction. It is not yet the solved state sensitivity; compare the linear solve in Section 10.
矩阵乘法也可以写成方向导数
A matrix product can be evaluated as a directional derivative
令 ,暂时固定设计输入 ()。向量 就是状态的切向量 :它规定所有未知量同时沿哪个方向变化。对每个残差关于 求导,就得到残差的切向量 。这只需要沿原来的计算过程应用链式法则,不必先存下整个矩阵;这就是 matrix-free(不显式组装矩阵)的含义,也就是自动微分的前向(tangent)模式:按原计算顺序,为每个中间量 算出它的切向量 。
Set while holding the design input fixed (). The vector is the state tangent : it specifies the direction in which all unknowns change together. Differentiating each residual with respect to gives the residual tangent . Apply the chain rule along the original computation without storing the full matrix: this is matrix-free evaluation, the forward (tangent) mode of automatic differentiation, which computes the tangent of every intermediate in the order of the original computation.
记号:上点 与上横线 是对离散程序求的前向与反向导数,与 Tapenade 等 AD 工具的变量 、 一一对应,下面的伪代码就这样命名。变分符号 只在附录 A 中出现,表示连续函数的变分。
Notation: the dot and the bar are forward- and reverse-mode derivatives of the discrete program. They correspond one to one to the variables and of AD tools such as Tapenade, and the pseudocode below uses these names. The variation symbol appears only in Appendix A, for variations of continuous functions.
先看最小的一步:。分别对两个输入求偏导,再乘各自的切向量并相加,就是 。在二阶格式中, 又来自重构,所以必须继续用链式法则传播到模板上的每个 。
Start with the smallest operation, . Differentiate with respect to each input, multiply by its tangent, then add to obtain . At second order, the trace tangents themselves come from reconstruction, so continue the chain rule to every stencil value .
手算一个面:取 ,则 、。给切向量 ,得到 。左残差的切向量加 ,右残差的切向量减 。这是一组局部演算数据,不是下方收敛流场中的固定面值。
Work one face by hand: with , obtain , . The tangents give . Add this to the tangent of the left residual and subtract it from the right one. These are illustrative local values, not fixed values of a face in the converged field below.
从状态切向量到通量切向量
From state tangents to flux tangents
令状态的切向量 ,设计固定()。无斜率限制器时,重构矩阵与状态无关,因此 、。链式法则必须经过重构、通量和分发三步。
Let the state tangent be with the design fixed (). Without a limiter, reconstruction matrices are state independent, so and . The chain rule must pass through reconstruction, flux evaluation and residual scatter.
二阶时,一个面依赖哪些未知量
Which unknowns a second-order face depends on
对一个内部二阶面,重构系数和该面通量的完整导数如下。这里 在重构面值上计算,不能在未经重构的中心值上计算。
For an interior second-order face, the reconstruction and complete flux derivative are as follows. Evaluate at the reconstructed face states, not at the unreconstructed centre values.
边界附近的二阶迹用单边斜率:第一内部面()左迹 ,最后内部面右迹 ;出口节点就在 ,,不需要重构。一阶时每个迹都取所在节点的值。出口物理通量的切向量为 。格点最后覆盖入口行,它的切向量是 。下面的播放器使用这些边界模板;上面的四点式只用于内部面。
Near the boundaries the second-order traces use one-sided slopes: the left trace of the first interior face () is and the right trace of the last interior face is ; the outlet node lies at , so needs no reconstruction. At first order each trace is the value of its node. The outlet physical-flux tangent is . For nodes, the inlet row is overwritten last, and its tangent is . The player below uses these boundary stencils; the four-point formula above applies only to interior faces.
播放器中的 是为了观察乘积而指定的测试种子 ,还不是第 10 节要求解的物理前向解。若选择单位向量 ,结果就是 Jacobian 的第 列;一般 是这些列的线性组合。伪代码就是 AD 工具生成的前向代码 的样子:输入种子 ,输出 。
The player uses the test seed with to illustrate a product. It is not the physical tangent solution of Section 10. Choosing would return column of the Jacobian; a general vector combines its columns. The pseudocode is what an AD tool generates as forward code : input seed , output .
逐句求导:残差代码的前向语句
The residual code differentiated statement by statement
前向模式对第 5 节残差代码的每一条赋值 用一次链式法则,得到 ,放在原语句之后、按同样的顺序执行;本节末尾的推导把这条规则用到面循环的四条语句上。下表把它用到本页的每一条语句,同时允许状态种子 与设计种子 :
Forward mode applies the chain rule once to every assignment of the residual code of Section 5, giving , placed after it and executed in the same order; the derivation at the end of this section applies the rule to the four statements of the face loop. The table applies it to every statement of this page's code, with a state seed and a design seed both allowed:
| # | 原始语句Primal statement | 前向语句(同一顺序)Tangent statement (same order) | 用到的局部导数Local derivative used |
|---|---|---|---|
| 1 | — | ||
| 2 | |||
| 3 | ;出口面 ; outlet face: | ;出口面 ; outlet face: | |
| 4 | (物理通量 ,) (physical flux: , ) | ||
| 5 | |||
| 6 | |||
| 7 | |||
| 8 |
前向代码为什么按原顺序执行,见本节末尾的推导;除了当前面上的原始值 ,不需要额外保存任何东西。设计只在两处进入:入口状态(格点页第 8 条约束)和源项(第 7 条)。取 、,这段代码返回 ;取 、,返回 的一列 ,也就是第 10 节切向方程 里的向量 (方程的右端是 ;这个 是向量,与面的下标 无关)。
Why the tangent code runs in the order of the primal code is explained in the derivation at the end of this section; nothing has to be stored beyond the primal values at the current face. The design enters in two places only: the inlet state (the constraint of statement 8 on the node page) and the sources (statement 7). With and the code returns ; with and it returns one column of , the vector of the tangent equation of Section 10 (whose right-hand side is ; this is a vector, not the face label ).
用本页的数值走一遍一个面:原始与前向
One face with the page’s own numbers: primal and tangent
取默认网格()、设计种子 ( 增加一个单位), 取第 10 节前向求解的收敛解,看面 (,左右是节点 4、5)。下表是它的原始量与切向量;二阶时每个迹由三个未知量重构,所以 也各由三个 组合。本节末尾的推导取节点 3、4 之间的面,设计种子是 ( 是入口值 ):
Take the default mesh (), the design seed (a unit increase of ) and for the converged solution of the tangent solve of Section 10, and look at face (, between nodes 4 and 5). The table lists its primal and tangent values; at second order each trace is reconstructed from three unknowns, so each combine three as well. The derivation at the end of this section works the face between nodes 3 and 4 with the design seed , where is the inlet value :
| 量Quantity | 公式Formula | 一阶First order | 二阶Second order |
|---|---|---|---|
所有面的 放在一起还是一个检查。 满足切向方程时,每个控制体的切向残差 都为零,所以面通量的切向量只在注入源项的控制体处跳变 ,其余处处相等。这里种子只加在节点 1上,于是一阶时 在第一个面上是 ,其后每个面都是 ,差值正好是 。图 13用中间的源项 把这个过程做成动画。
The of all faces together make one more check. When solves the tangent equation, the tangent residual of every control volume vanishes, so the face-flux tangent jumps by only where the seed injects source and is equal everywhere else. Here the seed acts on node 1 only, so at first order is on the first face and on every face after it: the difference is exactly . Figure 13 animates this for the middle source .
局部导数怎样组成全局 Jacobian
How local derivatives form the global Jacobian
固定网格和当前可微分支,把一个面的所有输入记为 ;它包括两侧的重构状态(格点页没有入口面:入口值 不进入任何面,只出现在被替换的入口行里)。令 、 分别是 对存储状态和设计变量的导数, 是局部通量 Jacobian。于是每个残差行的导数按同一个有向面关联矩阵 装配,这里的求和包括相应的边界通量映射(本页的这些求和是面循环的,即覆盖之前的残差 的);被替换的边界行另按第 9 节处理。前向模式按计算顺序传播切向量,状态种子是 ,设计种子是 :
At fixed mesh and differentiable branch, let collect all face inputs: the reconstructed states of both sides (the node page has no inlet face: the inlet value enters no face and appears only in the replaced inlet row). Let and be the derivatives of with respect to the stored states and the design variables, and the local flux Jacobian. Assemble each residual derivative with the same oriented face incidence , including the boundary flux maps (on this page these sums are those of the face loop, that is of the residual before any row is overwritten); replaced boundary rows are handled in Section 9. Forward mode carries the tangents in the order of the computation, with the state seed and the design seed :
无体源项时 。本页的设计变量不改变几何:源项进入 ;入口值不进入任何面,只出现在被替换的入口行里。二阶重构使 连接更多邻居,因此精确 Jacobian 的模板变宽;matrix-free 仍只需依次执行这些局部乘积,不必存储全局矩阵。
Set when there is no volume source. The design variables of this page do not change the geometry: the sources enter ; the inlet value enters no face and appears only in the replaced inlet row. Second-order reconstruction connects to more neighbours and widens the exact Jacobian stencil. Matrix-free evaluation performs these local products without storing the global matrix.
从程序到切向量:面循环的四条语句逐条线性化
From program to tangent: linearise the four statements of the face loop one by one
上面的 、 与 来自同一段短程序:残差代码对每个面依次执行四条语句。设 是面 读取的存储位置(两侧的 、,以及重构需要时它们的邻居),把这些位置上的状态排成 :
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 the reconstruction 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 face states and (the two arguments of this page's flux ; at first order ). (3) Flux: . (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 乘积就是按面执行这几次局部运算,不存储任何全局矩阵。下面的面板按这四条语句逐条播放本页一个面的数值。
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. The panel below plays the four statements with the numbers of one face of this page.
本页一个面的算例:前向扫描
A worked example on this page: the forward sweep of one face
取本页面循环里节点 3 与 4 之间的内部面 ();一阶重构只读两侧的状态,。设计种子取 , 是入口值 ;状态种子取本页前向求解对同一方向给出的解 ,于是下面每个切向量都是对 的导数。表中的数由本页代码在构建时算出:基准状态、前向求解与局部导数都用本页自己的程序。
Take the interior face of this page's face loop between nodes 3 and 4 (); first-order reconstruction reads only its two sides, . Seed the design with , where is the inlet value , 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 | 1.167 | 0.8688 | |
| 1.218 | 0.8349 | ||
| 设计种子design seed | 1 | 1 | |
| (2)重构(2) reconstruct | 1.167 | 0.8688 | |
| 1.218 | 0.8349 | ||
| (3)通量(3) flux | 0.6708 | 1.043 | |
| (4)分发(4) scatter | 0.6708 | 1.043 | |
| −0.6708 | −1.043 |
不直接进入这个面:, 全部来自状态的切向量。 不进入本页的任何面,;它只出现在被替换的入口行 里,所以 只有第 0 行非零。
does not enter this face directly: , and comes entirely from the state tangents. enters no face of this page, ; it appears only in the replaced inlet row , so has a single nonzero entry, in row 0.
所有面与源项的这些贡献加起来是替换之前的残差的切向量 。它在保留的行上于收敛解处为零;在被替换的行上它不为零(那一项是 ),因为残差代码在面循环之后把这一行整行覆盖,覆盖之后的切向量是约束的切向量。下面前向求解的方程是整个方程组的:
Summed over all faces and the source terms, these contributions give the tangent of the raw residual, . In the retained rows it vanishes at the converged solution; in the replaced row it does not (its entry there is ), because the residual code overwrites that row after the face loop and the tangent after the overwrite is the constraint's. The tangent equation solved further on is the system's:
按程序顺序执行四条语句,箭头上是各自的初等 Jacobian;下方是它们的乘积与形状(行 × 列,本页这个面的尺寸)。
The four statements in program order, each arrow labelled with its elementary Jacobian; below, their product with the shapes (rows × columns, this face of this page).
种子:( 是入口值 ),。在 上: 0.8688, 0.8349。
Seeds: ( is the inlet value ) and . On : 0.8688, 0.8349.
(1)取值:,只挑出 上的 2 个状态的切向量(上一步列出的数)。
(1) Gather: picks the tangents of the 2 states of (the numbers of the previous step).
(2)重构:。 0.8688, 0.8349。: 不进入这个面。
(2) Reconstruct: . 0.8688, 0.8349. : does not enter this face.
(3)通量: 1.043。
(3) Flux: 1.043.
(4)分发: 1.043,,即加上 −1.043。
(4) Scatter: 1.043 and , that is −1.043 added.
伴随:把同一个面循环倒过来
Adjoint: reverse the same face loop
本节实验:只给一个输出非零权重,预测反向影响能传到哪些输入;共享输入应累加所有路径。
Experiment: Seed one output and predict which inputs receive reverse contributions. Shared inputs accumulate every path.
先转置一个面的局部矩阵
First transpose the local matrix of one face
先忽略重构,一个内部面只有左右两个输入与两个残差贡献。前向局部矩阵如下:第二行是第一行的相反数,因为相同通量在右控制体里带负号。把这个小矩阵转置后作用在残差种子上,就自然得到先作差、再乘通量偏导的伴随步骤。
First omit reconstruction. One interior face has two inputs and two residual contributions. The second row of its forward matrix is the negative of the first because the right volume receives the negative flux. Applying the transpose of this small matrix to the residual seeds naturally gives the adjoint operations: subtract the seeds, then multiply by flux derivatives.
反向代码按与前向相反的顺序执行,每个输入的横线量都用 += 累加:由种子 、 先得到 、,再得到 、。给定第 8 节的局部数值,再取 ,则 、。检查对偶恒等式的两边:,而 。这就是转置关系在一个面上的具体含义。
The reverse code runs in the opposite order of the forward code and accumulates the bar of every input with +=: from the seeds and it forms and , then and . Use the local values from Section 8 and : then , and . Both sides of the duality identity agree: and . This is the transpose identity on a single face.
为什么二阶必须累加到远邻
Why second order reaches beyond the two neighbours
如果 ,反向传播就给这三个未知量的伴随分别累加:、、。右迹也做同样的事。某个未知量可能从多个面、同一面的两条重构路径收到贡献,因此必须使用“+=”,不能后到的贡献覆盖先到的贡献。
For , reverse propagation adds to the adjoints of those three unknowns: , and . Apply the same rule to the right trace. An unknown can receive contributions from several faces and from both reconstruction paths on one face, so use += rather than overwriting earlier contributions.
从残差权重返回状态权重
From residual weights back to state weights
伴随是同一条计算链按相反顺序的转置。以残差种子 出发:先收集每个面上通量的伴随 ;再乘通量偏导,得到左右迹的伴随 、;最后通过 累加(+=)到所有受影响未知量的伴随 。不能只把结果放回面两侧的中心。
The adjoint is the transpose of the same chain, run in reverse order. Start from the residual seed : gather the adjoint of each face flux, ; multiply by the two flux derivatives to get the trace adjoints and ; then accumulate (+=) through into the adjoint of every affected unknown. Returning contributions only to the two adjacent centres would omit reconstruction derivatives.
图 11把 与 写成矩阵形式对照。前向模式一次扫描得到 的各列按 加权之和:种子 恰好取出第 列 。反向模式一次扫描得到各行按 加权之和:种子 取出第 行 。所以前向模式适合“少输入、多输出”,反向模式适合“多输入、少输出”——一个目标对很多设计变量求导,正是后者。两种模式都不组装矩阵;图中的矩阵由本页的 Jacobian 画出,只为显示每一步用到哪些元素。
Figure 11 sets beside in matrix form. One forward sweep delivers the columns of weighted by : the seed picks out exactly column , . One reverse sweep delivers the rows weighted by : the seed picks out row , . Forward mode therefore suits few inputs and many outputs, reverse mode many inputs and few outputs, which is the case of one objective and many design variables. Neither mode assembles the matrix; the figure draws it from this page's Jacobian only to show which entries each step uses.
入口约束怎样反向求导
Differentiating the inlet constraint in reverse
格点强边界按相反顺序处理:前向代码最后才用约束 覆盖第 0 行,所以严格按反序它的反向语句最先执行:用掉种子 (,因为 ),并令被替换行的种子 (面循环写入该行的通量已被覆盖,不再得到权重)。 也可以放在完成所有面及重构的反向累加之后:+= 可以交换次序,只要 在面循环读取它之前已经清零;本节的播放器就这样做。这一步必须累加,不得覆盖已经从其他残差行传回 的贡献。
For a strong node-based boundary, reverse the order: the forward code overwrites row 0 last with the constraint , so in strict reverse order its statement runs first: it uses the seed (, because ) and sets the replaced row's seed (the face-loop flux written into that row was overwritten and receives no weight). The statement may equally come after all face and reconstruction contributions: += commutes, as long as is cleared before the face loop reads it, and the players of this section do so. This step must accumulate; do not overwrite contributions already propagated to from other residual rows.
点积检验 ,即 AD 的对偶恒等式 ,检查转置关系,但两个方向若同时漏掉重构,仍可能通过。第 14 节会用独立复数残差和总导数再做检验。
The dot test , the AD duality identity , checks transposition, but can pass if both operators omit reconstruction. Section 14 additionally checks an independent complex residual and total derivatives.
“倒过来”指依赖链的次序反转,不要求把相互独立的面编号倒序遍历。这里每个面读取固定的 与种子 ,并向独立的输出数组 累加,正序或倒序遍历面都给出同一线性算子(浮点舍入次序可有微小差异)。GS 的行更新会读刚写入的新值,因而其反向扫描是另一件必须明确处理的事。伪代码就是 AD 工具生成的反向代码 的样子:种子 ,输出 。
Reversal concerns dependency order, not necessarily the numbering of independent faces. Here each face reads the fixed and seed and adds to a separate output array ; either face order gives the same linear operator up to rounding. GS row updates read freshly written values, so reversing their scan is a separate issue that must be treated explicitly. The pseudocode is what an AD tool generates as reverse code : seed , output .
逐句倒推:残差代码的反向语句
The residual code reversed statement by statement
下表从最后一条原始语句往回,列出第 8 节表中每条语句的伴随语句 ,种子是残差的伴随 ;面循环里,每个面内的语句也倒序执行。为什么倒序、为什么用 +=,见本节末尾的推导:
The table lists the adjoint statement of every statement of the table in Section 8, from the last primal statement back to the first, with the residual adjoint as the seed; inside the face loop the statements of each face also run in reverse. Why in this order and why with +=: see the derivation at the end of this section:
| # | 原始语句Primal statement | 伴随语句(倒序执行)Adjoint statement (run in reverse order) |
|---|---|---|
| 8 | ||
| 7 | ||
| 6 | ||
| 5 | ||
| 4 | ||
| 3 | ;出口面 ; outlet face: | ;出口面 ; outlet face: |
| 2 | ||
| 1 | 无: 是种子nothing: is the seed |
汇总的是加权输出 对 的敏感度。为什么倒序、为什么每个横线量都用 += 累加,本节末尾的推导按面循环的四条语句给出(程序里的一次扇出,在反向扫描里变成扇入)。本页的计数:一阶时 被两个迹读取,二阶时被四个面上的六个迹读取;格点页的 被两个控制体的源项读取。覆盖赋值则切断旧值的所有路径:第 8 条把第 0 行整个改写,所以它的伴随语句用完 后令 ,此前写入第 0 行的通量和源项不再得到权重。本节的播放器把 放在面循环之后执行;+= 可以交换次序,结果相同,关键是面循环读取 之前它已经清零。
collects the sensitivity of the weighted output to . Why the statements run in reverse order and why every bar accumulates with += is derived for the four statements of the face loop in the derivation at the end of this section (a fan-out of the program becomes a fan-in of the reverse sweep). The counts on this page: at first order is read by two traces, at second order by six traces on four faces, and on the node page by the sources of two volumes. An assignment that overwrites a variable cuts every path through its old value: statement 8 rewrites row 0 completely, so its adjoint statement uses and then sets ; the flux and the source written into row 0 before receive no weight. The player of this section adds into after the face loop instead; += commutes, so the result is the same, as long as is cleared before the face loop reads it.
每条语句都满足一个局部的对偶关系:对 ,有 。把所有语句串起来,就是整段代码的对偶恒等式( 时就是上面的点积检验),本节末尾的推导对面循环逐级证明它:
Every statement satisfies a local duality: for , . Chained over all the statements this is the duality identity of the whole code (with it is the dot test above), which the derivation at the end of this section proves step by step for the face loop:
同一个面倒过来:伴随值与对偶检查
The same face in reverse: adjoint values and the duality check
继续第 8 节的面 ,种子取第 11 节伴随方程的解 。反向语句 6、5 先把两侧的种子收集成 ,语句 4 再乘上同一对偏导:
Continue with the face of Section 8, now seeded with the solution of the adjoint equation of Section 11, . The reverse statements 6 and 5 first gather the seeds of the two sides into , then statement 4 multiplies by the same pair of partial derivatives:
| 量Quantity | 公式Formula | 一阶First order | 二阶Second order |
|---|---|---|---|
| 输出端output side | |||
| 输入端input side |
两端一致到舍入误差:局部对偶关系在这个面上成立。它对整段代码同样成立,而且此时的数值有明确含义: 满足切向方程,所以 ; 满足伴随方程,所以 。代入对偶恒等式得 ,即 :伴随给出的梯度分量等于前向给出的方向导数。一阶时两者都是 (二阶 )。本节末尾的推导对节点 3、4 之间的面取同样的种子 与设计种子 (),逐级列出内积。
The two sides agree to round-off: the local duality holds on this face. It holds for the whole code as well, and there the numbers mean something: solves the tangent equation, so ; solves the adjoint equation, so . The duality identity then reads , i.e. : the gradient component from the adjoint equals the directional derivative from the tangent. At first order both are (at second order ). The derivation at the end of this section reverses the face between nodes 3 and 4 with the same seed and the design seed (), listing the inner products stage by stage.
图 12把点积检验拆成三层:残差层 、面层 (格点页再加约束行的 ) 与状态层 ,种子与本节播放器相同,、。面层等于残差层:格点页的第 0 行被约束 替换,面循环写入它的量被丢弃,约束行自己的切向量是 ,所以 ,面上的 用清零了第 0 行的种子 ( 取 0)。面层因此是各面之和 加上约束行的 ,共 ,与残差层相等;状态层等于面层,因为每个面上 ,而 正是反向语句 在状态层里的一项。逐项看,三层的各项并不相同;逐层求和,三个总和相等。
Figure 12 splits the dot test into three levels: the residual level , the face level (on the node page plus the constraint row's ) and the state level , with the seeds of this section's players, and . The face level equals the residual level: row 0 of the node page is replaced by the constraint , so what the face loop wrote into it is discarded and the constraint row's own tangent is ; thus , and the face terms use with the seed that has row 0 cleared ( taken as 0). The face level is therefore the sum over the faces, , plus the constraint row's , in all, equal to the residual level. The state level equals the face level because on every face, and is the term that the reverse statement contributes to it. Term by term the three levels differ; summed level by level, the three totals are equal.
为什么反向累加给出矩阵转置
Why reverse accumulation gives the matrix transpose
把所有残差和状态分量展平后,下面只是交换两个有限求和的次序。残差种子 是任意给定的,尚未求解伴随方程;只有选择满足第 11 节线性方程的 ,才得到该目标对应的离散伴随解。取 、,前向得到 ;以 为种子反向扫描,得到 。下式就是 AD 的对偶(点积)恒等式:
After flattening all residual and state components, the following simply exchanges two finite sums. The residual seed is arbitrary and has not solved an adjoint equation; only satisfying Section 11 gives the discrete adjoint solution for the chosen objective. With and the forward sweep gives ; the reverse sweep seeded with gives . The identity below is the duality (dot-product) identity of AD:
格点页的残差代码在面循环与源项之后还有一条语句(5):把入口节点的守恒行覆盖为约束 。 是保留行的对角掩码, 是第 0 个单位向量;覆盖之前的残差、Jacobian 与设计 Jacobian 记作 、、,面循环与源项写出的就是它们。
The residual code of the node page has one more statement (5) after the face loop and the sources: it overwrites the balance row of the inlet node by the constraint . masks the retained rows and is the first unit vector; the residual, Jacobian and design Jacobian before the overwrite are , and , and the face loop and the sources write exactly these.
这条语句最后执行,所以它的反向语句(5′)最先执行:种子 分成两部分,保留行的部分 交给转置的面循环与源项,第 0 行的部分 交给约束,令 (因为 )与 (因为 )。完整转置必须包含约束的 Jacobian,而不是简单地把入口伴随设为零;被替换的第 0 行也不能把 交给面循环。
This statement runs last, so its reverse statement (5′) runs first: the seed splits into the part of the retained rows, , which goes to the transposed face loop and the sources, and the part of row 0, , which goes to the constraint: (because ) and (because ). The transpose includes the constraint Jacobian; it does not simply zero the inlet adjoint, and the replaced row 0 must not hand to the face loop.
内部面通量向左加、向右减,所以先把这两条路径合并成通量的伴随 ;之后按通量、重构、原始状态的相反依赖顺序传播。反向语句中的 += 很重要:同一状态会被多个面和多个求积点读取,这些贡献必须累加。这里没有分部积分,也没有另行猜测连续伴随边界条件。
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 +=:
从切向量到伴随:转置把语句顺序反过来
From tangent to adjoint: the transpose runs the statements backwards
一个面的切向量是四个初等 Jacobian 的乘积。乘积的转置等于各个因子的转置按相反的顺序相乘:
The tangent of one face is the product of four elementary Jacobians. The transpose of a product is the product of the transposes in the reverse order:
仍从右往左读,这就是反向扫描的语句顺序:程序最后执行的分发最先处理,最先执行的取值最后处理。先把所有伴随置零,再以 为种子逐条执行;覆盖(5)在整个面循环之后执行,所以它的反向语句(5′)排在最前:
Read again from right to left, this is the order of the reverse sweep: the scatter, which the program ran last, comes first, and the gather, which it ran first, comes last. Zero every adjoint, seed with , and run; the overwrite (5) ran after the whole face loop, so its reverse statement (5′) comes first:
面循环结束后,源项的一对也按同样的规则反过来:,,;被替换的行已由上面的(5′)处理:它们的种子不进入面循环,写进这些行的源项因此得不到权重。本页代码里每一条语句的伴随语句(顺序与前向相反)也都在本节里逐句写出。
After the face loop the source is reversed by the same rule: , , ; the replaced rows were handled by (5′) above: their seeds do not enter the face loop, so a source written into such a row receives no weight. The adjoint statement of every statement of this page's code, in the reverse order, is written out in this section as well.
分发与取值互为转置。分发 把一个面通量写进两行;它的转置 从这两行读回一个数,是一次取值。取值 从全局数组读出几个状态;它的转置 把它们的伴随加回原来的位置,是一次带累加的分发。
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.
对偶恒等式的逐步证明,以及它的数值形式:点积测试
The duality identity step by step, and its numerical form: the dot-product test
把同一个面的前向语句与反向语句配对。每一步只用一次转置的定义 ,把一个初等 Jacobian 从内积的右边移到左边,所以内积在每一级都相同:
Pair the forward and reverse statements of the same face. Each step uses the definition of the transpose once, , to move one elementary Jacobian from the right of the inner product to the left, so the inner product is the same at every stage:
对所有面求和,再加上源项的一对 (),右边的两项正好累加成反向扫描的输出 与 :
Summed over all faces, with the pair of the source added ( with ), the two terms on the right accumulate exactly into the outputs and of the reverse sweep:
被替换的行再加上它们约束的一对 ( 是切向量一节里覆盖语句的切向量):它正是(5′)里约束的状态路径与设计路径对 的贡献;所以同一个恒等式对本页的方程组也成立( 与 取整个反向扫描结束时的值),下面的测试用的就是它:
The replaced rows add the pair of their constraints, ( is the tangent of the overwrite in the tangent section): it is exactly what the state path and the design path of the constraints in (5′) contribute to . The same identity therefore holds for this page's system ( and taken at the end of the whole reverse sweep), and the test below uses that system:
这就是 AD 的点积测试:任取种子 与 ,前向扫描算出 ,反向扫描算出 ,两边之差应在舍入误差量级;只要有一条反向语句的转置写错,差就是 的。本页在构建时用它自己的 matrix-free 前向与反向乘积( 个分量的确定性种子,设计方向取下面算例的 )做了这个测试:,,相对差 。
This is the dot-product test of AD: take any seeds and , run the forward sweep for and the reverse sweep for ; the two sides must agree to round-off, and a single wrong transpose in the reverse statements gives an difference. This page ran the test on its own matrix-free forward and reverse products at build time (deterministic seeds with components, the design direction the of the example below): and , a relative difference of .
同一个面的反向扫描
The reverse sweep of the same face
反向种子取本页伴随求解的解 ;面循环的种子是掩码后的 ,这个面的两行都是保留行,所以它在这两行等于 。于是 与 是这个面对 与 的贡献(也是对 与 的贡献)。下表按反向语句的顺序列出伴随值;最后一列是每一级的内积:先是 ,然后是 、,最后是 ,切向量取第 8 节“本页一个面的算例:前向扫描”里同一个面的前向算例的值。
Seed the reverse sweep with this page's adjoint solution, ; the seed of the face loop is the masked one, , and both rows of this face are retained, so it equals in them. Then and are this face's contributions to and (and so to and ). The table lists the adjoints in the order of the reverse statements; the last column is the inner product at each stage: first , then , , and last , with the tangents of the forward example of the same face, the one under “A worked example on this page: the forward sweep of one face” in Section 8.
| 语句Statement | 量Quantity | 伴随Adjoint | 内积Inner product |
|---|---|---|---|
| (0)置零(0) zero | 0 | ||
| (4′)种子(4′) seed | 0.05619 | 0.005043835575 | |
| 0.05135 | |||
| (4′)取值(4′) gather | 0.004837 | 0.005043835575 | |
| (3′)通量的转置(3′) flux, transposed | 0.006693 | 0.005043835575 | |
| −9.233×10⁻⁴ | |||
| (2′)(1′)加回(2′) (1′) add back | 0.006693 | 0.005043835575 | |
| −9.233×10⁻⁴ | |||
| 0 |
4 级内积相同,最大相对差 :对偶恒等式在这些数上成立;下面的面板按反向语句逐条播放这些数。这个面不直接依赖 ,所以 ; 的伴随来自它真正进入的那些行。
The 4 inner products agree, with a largest relative difference of : the duality identity holds in these numbers. The panel below 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.05619, 0.05135(这个面的两行都是保留行,掩码后的种子 在这两行等于 )。
(0) Zero and seed: , ; 0.05619, 0.05135 (both rows of this face are retained, so the masked seed equals in them).
(4′)分发的转置是一次取值: 0.004837。内积 0.005043835575。
(4′) The transposed scatter is a gather: 0.004837. Inner product 0.005043835575.
(3′)通量的转置:, 0.006693, −9.233×10⁻⁴。内积 0.005043835575。
(3′) The transposed flux: , 0.006693, −9.233×10⁻⁴. Inner product 0.005043835575.
(2′)重构的转置:,,这个面给出 0()。
(2′) The transposed reconstruction: and ; this face gives 0 ().
(1′)取值的转置是带累加的分发:, 0.006693, −9.233×10⁻⁴。内积 0.005043835575。
(1′) The transposed gather is a scatter with accumulation: , 0.006693, −9.233×10⁻⁴. Inner product 0.005043835575.
对偶:4 级内积相同(最大相对差 0), 在这个面上逐级成立。
Duality: the 4 inner products agree (largest relative difference 0); holds stage by stage on this face.
前向求解:设计改变后,解怎样变化
Forward solve: how the solution changes with the design
从方向乘积到一个需要求解的方程
From a directional product to an equation to solve
前面可以对任意向量 计算 。现在给设计输入一个种子(方向):状态也随之变化,它的切向量是 ,且新状态仍要满足 。对这个恒等式前向求导(多元链式法则),得到两部分:经状态切向量的间接影响与设计直接进入残差的影响。
We can now evaluate for any supplied . Now seed the design input with a direction : the state changes too, with tangent , and must still satisfy . Differentiating this identity in forward mode (the multivariable chain rule) gives two terms: the indirect effect through the state tangent and the direct design effect in the residual.
例如取种子 ,即第 个源区间增加一个单位的方向:格点页上它进入相邻两个半控制体,直接项 ( 时只有 :第 0 行是入口约束,不含源项),其余分量为零;随后 的多个分量会响应。这区分了“输入在哪里施加”和“影响传播到哪里”。先固定收敛的原始解,再求这个线性系统,最后用 得到目标的方向导数。
For example, the seed , a unit change in source interval , enters the two adjacent half volumes on the node page, giving direct terms (for only : row 0 is the inlet constraint and carries no source) with other components zero, yet many state components can respond in . This separates where the input is applied from where its effects propagate. Freeze the converged primal state, solve the linear system, then obtain the objective directional derivative as .
沿设计种子 对收敛残差 前向求导,令 ,得到精确 tangent 方程,它的解是状态切向量 。在该线性求解期间固定原始解、 和一阶预处理矩阵 ;修正量 是迭代的增量,不是导数。
Differentiate the converged residual forward along the design seed and set . The exact tangent equation below has the state tangent as its solution. The primal state, exact and first-order preconditioning matrix remain fixed during this linear solve; the correction is an iteration increment, not a derivative.
用一阶 Jacobian 求二阶方程的修正
A first-order Jacobian for corrections to the second-order equation
每次修正方程从零开始,用所选次数的 GS 扫描,得到线性算子 。一次扫描时 。精确求出修正方程时,可以把这一缺陷校正写成 ();有限次 GS 时应使用下面的实际迭代式。
Starting each correction from zero and taking the selected fixed number of GS sweeps defines a linear operator . One sweep gives . If the correction is solved exactly, this defect correction becomes for ; with finite GS sweeps, the actual iteration is the following.
收敛需要误差传播矩阵的谱半径小于 1;GS 或低阶预处理并不对任意二阶问题保证收敛。下面可更换设计变量,查看解及真实线性残差。
Convergence requires a spectral radius below 1 for the error propagator; GS or low-order preconditioning does not guarantee convergence for arbitrary second-order problems. Select a design variable below and inspect the solution and true linear residual.
前向固定点求解:每次检查精确线性残差
Forward fixed point: check the exact linear residual each time
f = design_derivative(u, D_d) # D_d is the design seed Ḋ
v = zeros(M); P = first_order_jacobian(u)
repeat:
r = forward_product(u, v) + f # A v + f, one forward sweep
if norm(r) <= tolerance: return v
delta = GS(P, -r, zero_start, fixed_sweep_count) # the correction Δv
v = v + omega*delta一次 GS 应用把右端向量变成修正,记为线性映射 。将更新式代入误差 ,就得到误差传播矩阵;其谱半径小于 1 才保证对任意初始误差收敛。残差曲线可以观察本次求解,但不能替代一般收敛条件。
One GS application maps a right-hand side to a correction; call this linear map . Substitute the update into to obtain the error propagation matrix. A spectral radius below 1 guarantees convergence for every initial error. A residual plot describes this solve, not a general convergence theorem.
残差对设计变量的导数:逐列写出
The residual’s design derivative, column by column
切向方程 里的向量 (右端是 )需要每一行残差在状态固定时对每个设计分量的偏导。读取设计的只有两类语句:源项积分,以及入口值(格点页是约束行)。所以 的每一列都很稀疏:
The vector of the tangent equation (its right-hand side is ) needs the partial derivative of every residual row with respect to every design component, with the state held fixed. Only two kinds of statement read the design: the source integrals, and the inlet value (the constraint row on the node page). So every column of is sparse:
在前向代码里,这两类列就是第 8 节表中取 的第 7 条语句()和第 8 条语句():设计种子直接改变的只有它自己的两个半控制体或约束行。
In the tangent code these columns are statements 7 () and 8 () of the table in Section 8 with : a design seed changes directly only its own two half volumes, or the constraint row.
图 13跟着一个种子走完代码与求解。取区间右端在 的源项 :先用 执行面循环,所有面的 都是零,只有源项语句让节点 3、4(各 )的 不为零,这就是 。接着,前向 GS 从第 0 行扫到最后一行:右端在种子上游为零,而下三角求解只把已经更新的上游值传给下一行,所以上游各行保持为零,从第 3 行起向下游逐行点亮——这正是信息沿特征速度 的传播方向。求解收敛后再以 执行一次面循环:每个面的 都亮起来,只在种子处跳变,最后所有 回到零。
Figure 13 follows one seed through the code and the solve. Take the source of the interval ending at . First the face loop runs with : every is zero, and only the source statement makes of nodes 3 and 4 ( each) nonzero; that is . Then the forward GS sweep runs from row 0 to the last row: the right-hand side is zero upstream of the seed, and the lower-triangular solve passes only already updated upstream values to the next row, so the rows upstream stay zero and from row 3 on they light up one after another downstream, the direction in which information travels with the characteristic speed . When the solve has converged, the face loop runs once more with : every face flux tangent lights up, jumps only at the seed, and at the end every is back at zero.
把隐函数求导展开成可求解的系统
Expand implicit differentiation into a solvable system
设原始残差已经收敛,且当前光滑分支上的 非奇异。取设计种子(方向),令 ;前向求解得到的状态切向量是 。对每一个残差分量逐项求导,零残差的切向量仍为零:
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.
实现分三步:状态种子 、设计种子 ,前向求值残差得到 ;用精确 的乘积求解 ;再以 和 为种子对目标前向求导,得到 。一次 只是算子调用(一次前向扫描),不是求解 。稳态导数的对象是收敛残差方程,不是 GS/GMRES 的迭代次数。在不可微点(如 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 GS/GMRES iteration count. At nondifferentiable points (a limiter switch or a branch test) these ordinary derivatives need separate treatment.
用面循环的语句来看:右端 是 时的一次前向扫描。它在保留的行上是面循环的 ,只有 与源项两条路径起作用;在被替换的行上是约束自己的设计路径,即第 0 行的 ,任何面循环里都没有它:对入口值 ,整个右端都来自这条路径(没有任何面读到它,)。迭代中的每个 是 时的一次前向扫描:保留的行是面循环的,被替换的行是约束对状态的路径。一次前向求解的全部工作就是这段面循环与覆盖以不同的种子反复执行。对 个设计分量逐个求导,要做 次这样的求解。
In terms of the face loop's statements: the right-hand side is one forward sweep with . In the retained rows it is the face loop's , in which only the paths and the source are active; in the replaced rows it is the constraint's own design path, in row 0, which no face loop contains: for the inlet value the whole right-hand side comes from that path (no face reads it, ). Every of the iteration is one forward sweep with : the face loop in the retained rows and the constraint's state path in the replaced rows. All the work of a tangent solve is this face loop and overwrite run again and again with different seeds, and the derivatives with respect to all design components take such solves.
伴随求解:从目标函数向回传递
Adjoint solve: propagate back from the objective
为什么转置方程能省去很多前向求解
Why a transpose solve avoids many forward solves
如果有很多设计变量而只关心一个 ,逐个计算 要解很多线性系统。利用 ,找一个 使 ,便可把 改写为只含 的内积。整个推导只用矩阵转置恒等式。
With many design variables but one objective, solving for each requires many linear systems. Given , choose satisfying . A transpose identity rewrites as an inner product with , eliminating .
这里的两个负号相消。也可以采用相反号的伴随变量,但方程和梯度公式必须一起改变。本页统一用 ,所以右端是 。伴随不是原始流场倒着求一次,而是一个由目标梯度驱动的转置线性问题。
The two minus signs cancel. An opposite adjoint sign convention is possible, but the equation and gradient formula must both change. This page uses , giving right-hand side . The adjoint is a transpose linear problem driven by the objective gradient, not the primal flow solved backwards.
图 14把代价直接数出来。用页首的线性求解器(默认 GS 固定点,每次一轮扫描),前向对每个设计分量解一次 ,每次只得到梯度的一个分量 ;伴随只解一次 ,再做一次乘积 就得到全部分量。默认设置下 9 次前向求解共用 210 次 GS 扫描,伴随只用 24 次;两条路线给出的 9 个分量相差不超过 。为什么伴随的代价与设计分量数无关,见本节末尾的推导与两条路线的代价表。
Figure 14 counts the cost. With the linear solver chosen at the top (by default the GS fixed point with one sweep per application), the tangent solves once per design component and obtains one gradient component each time; the adjoint solves once, and one product then gives every component. With the default settings the 9 tangent solves take 210 GS sweeps in total and the adjoint only 24; the 9 components of the two routes differ by at most . Why the adjoint's cost does not depend on the number of design components is derived, with the cost table of the two routes, at the end of this section.
ψ 的含义:给一个控制体多加一点源项
What ψ means: add a little source to one control volume
把控制体 的积分源项增加一个数 ,第 行变成 ,其余行不变。这等于求解 。按前向模式,给源项一个单位种子 (其余种子为零):残差的切向量 必须为零,所以 ;再用伴随方程 ,即 :
Increase the integrated source of control volume by a number : row becomes and the other rows are unchanged. This is the same as solving . In forward mode, give the source the unit seed (all other seeds zero): the residual tangent must vanish, so ; the adjoint equation , i.e. , then gives
所以 是“控制体 多注入单位积分源项时目标函数减少多少”。负号来自残差约定“流出减流入减源项”。若源项密度在宽度为 的控制体上均匀增加、种子为 ,则 ,于是 。格点页第 0 行是约束 :在第 0 行加 等于把入口值改为 ,所以 ,与第 12 节一致;源项不进入第 0 行。图 15不用伴随,直接测量这个关系:对每个 解两次完整的非线性方程 (步长 是一个数),用中心差分得到 。两者的最大相对差是 ,正是中心差分 截断误差的量级。
So is the decrease of the objective when one unit of integrated source is injected into control volume . The minus sign comes from the residual convention "outflow minus inflow minus source". If the source density rises evenly over a volume of width with the seed , then and . On the node page row 0 is the constraint : adding to row 0 changes the inlet value to , so as in Section 12; sources do not enter row 0. Figure 15 measures this relation without the adjoint: for every it solves the full nonlinear equations twice (the step is a number) and forms the central difference . The largest relative difference is , the size of the truncation error of a central difference.
图中的两个箭头说明方向。在 注入的源项沿流动改变下游的 ; 在每个控制体测量失配 (橙色柱),所以 汇总的是 下游的失配。出口附近下游几乎没有控制体, 很小;向入口方向每经过一个失配为负的控制体, 就增大,经过失配为正的控制体(这里在入口附近)则减小。这就是伴随“从目标函数向回传递”的含义,也是连续伴随从 出发、向入口积分的离散版本。
The two arrows give the directions. A source injected at changes downstream along the flow; measures the mismatch in every volume (orange bars), so collects the mismatch downstream of . Near the outlet almost nothing lies downstream and is small; moving toward the inlet, grows across every volume with a negative mismatch and shrinks across a volume with a positive one (here near the inlet). This is what "propagating back from the objective" means, and it is the discrete counterpart of integrating the continuous adjoint from toward the inlet.
采用拉格朗日函数 ,令 中状态切向量 的系数为零,得到 。修正方程使用 ,并按相反方向扫描;修正量记作 。
With , setting the coefficient of the state tangent in to zero gives . The correction uses , scans in the opposite direction and is written .
相同固定次数、零初值的反向扫描恰好构成 ,不是把矩阵转置后仍沿原方向扫描。于是实际更新和误差传播如下。
With the same fixed sweep count and zero initial correction, reverse sweeps implement . Transposing the matrix while retaining the original sweep direction would give a different preconditioner. The actual update and error propagation are therefore:
注意 并不直接等于 tangent 迭代矩阵的转置,因为乘法顺序不同。这里解的是残差伴随变量 ;若逆向微分整个迭代程序,需要额外说明其迭代伴随变量与 的转换。
Note that is not directly the transpose of the tangent iteration matrix: the multiplication order differs. We solve for the residual adjoint . Reverse differentiation of the entire iteration program would require an additional relation between its iterative adjoint variable and .
伴随固定点求解与反向 GS
Adjoint fixed point and reverse GS
g = objective_state_gradient(u)
psi = zeros(M); PT = transpose(first_order_jacobian(u))
repeat:
r = adjoint_product(u, psi) + g # A^T psi + g, one reverse sweep seeded with psi
if norm(r) <= tolerance: return psi
delta = zeros(M) # the correction Δψ
repeat fixed_number_of_sweeps:
for i = M-1,...,0:
delta[i] = (-r[i] - Σ(j≠i) PT[i,j]*delta[j]) / PT[i,i]
psi = psi + omega*delta一轮正向 GS 解下三角系统 ,所以 ;转置后是 ,这是上三角系统,需要从最后一行往前解。对固定的多轮扫描,同样需要转置整个零初值映射,而不只是矩阵名加一个 T。下方检验也比较 与 。
One forward GS sweep solves , so . Its transpose is , an upper-triangular solve performed from the last row backwards. With several fixed sweeps, transpose the entire zero-start map, not merely the matrix name. The page also checks .
把伴随方程逐行写出来
The adjoint equation written out row by row
转置把第 6 节的每个系数移到镜像位置。一阶时, 的第 行收集的是 的第 列的三个元素:
Transposition moves every coefficient of Section 6 to its mirror position. At first order, row of collects the three entries of column of :
按面重新分组,就是反向面循环算出的结果:第 行等于右面的 加上左面的 ,其中 是 跨过该面的差:
Regrouped by faces, this is exactly what the reverse face loop computes: row is the of its right face plus the of its left face, where is the jump of across that face:
一阶时这两个系数在同一个状态上计算:,,所以对角元是 。代入并除以 (),第 行成为
At first order the two coefficients are evaluated at the same state, and , so the diagonal entry is . Substituting and dividing by (with ), row becomes
令 :第一项是 的中心差分,第二项是 。所以离散伴随方程就是连续伴随方程 (第 16 节、附录 A)加上一个 的数值耗散,一阶格式的伴随因此在远离入口处一阶收敛到 。入口处不然:第 16 节的入口伴随边界层使入口第一个控制体的误差不随加密消失,全域 误差只有半阶。这一行还说明了方向:把它解成 ,下游邻居 的权重 在本页各内部行上介于 与 之间(入口一侧的第 1 行取 ,权重是 ): 主要由下游决定,所以从出口一侧开始的反向 GS 扫描收敛得快。
With the first term is a centred difference of and the second is . So the discrete adjoint equation is the continuous adjoint equation (Section 16, Appendix A) plus a numerical dissipation of order , and away from the inlet the adjoint of the first-order scheme converges to at first order. Not at the inlet: the inlet adjoint boundary layer of Section 16 keeps the error of the first control volume from vanishing under refinement, and the global error is only of order one half. The row also shows the direction: solved as , the weight of the downstream neighbour lies between and on the interior rows of this page (row 1, on the inlet side, takes and has the weight ): is set mainly by what lies downstream, which is why the reverse GS sweep that starts on the outlet side converges fast.
两端的行不同(下式)。约束行决定约束乘子 :它由右侧的伴随给出,就是第 12 节的 。第 1 行反而没有 的项:约束行在转置面循环之前已令 ,第一个面只传回 ,所以第 1 行是内部公式取 ;和格心页的入口单元一样,入口外的伴随被当作零,这也是入口伴随边界层的来源(第 16 节)。出口节点的行只含物理通量的 ,所以 ,对应连续伴随的终端条件 。
The end rows differ (below). The constraint row sets the constraint multiplier from the adjoint on its right; it is the of Section 12. Row 1, in turn, has no term in : the constraint row clears before the transposed face loop, so the first face passes back only , and row 1 is the interior formula with . As in the inlet cell of the cell page, the adjoint outside the inlet is taken to be zero, which is the source of the inlet adjoint boundary layer (Section 16). The row of the outlet node holds only the of the physical flux, so : the counterpart of the terminal condition of the continuous adjoint.
下面的动画用伴随 GS 固定点求解 (不论页首选择哪种线性求解器,这里都是 GS)。先逐行播放第一次修正:反向扫描从最后一行走到第 0 行, 从出口一侧开始出现,逐步延伸到入口。这与图 15的方向一致,也与连续伴随从 出发、向入口积分的方向相同。之后每帧是一次外层迭代,下方的误差柱是出口一侧先变短;与上面的原始动画一样,迭代越过停止准则,一直做到残差停在舍入平台上。
The animation below solves with the adjoint GS fixed point (always GS here, whichever linear solver is selected at the top). The first correction is played row by row: the reverse sweep runs from the last row to row 0, so appears on the outlet side first and extends toward the inlet. This is the direction of Figure 15, and the direction in which the continuous adjoint is integrated from toward the inlet. After that, each frame is one outer iteration, and the error bars below shrink on the outlet side first; as in the primal animation, the iteration continues past the stopping criterion until the residual rests on its round-off floor.
从有限维 Lagrangian 逐步消去状态敏感度
Eliminate state sensitivity from a finite-dimensional Lagrangian
在满足 的解上 。先把 固定,沿切向量 、 对 前向求导,再按 和 收集系数。即使 随设计改变,额外的 项也在收敛解上为零。
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 .
精确算子决定求解哪个方程;GS 固定点、GS 预处理 GMRES 与 LU(附录 C)决定如何求解。较便宜的近似 Jacobian(例如一阶或冻结部分系数的 Jacobian)可以用于迭代修正或预处理,但不能替代精确转置。应监控真残差 。若其不为零,当前伴随与精确 tangent 的方向导数差为 ,因此只看迭代次数或预处理残差不够。
The exact operator determines the equation; the GS fixed point, GS-preconditioned GMRES and LU (Appendix C) determine 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.
为什么一次伴随求解能代替 次前向求解
Why one adjoint solve replaces tangent solves
同样的结论也可以不经过 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) | 次 |
本页 、:前向路线要解 9 个 9 阶的线性方程组,伴随路线只解 1 个。设计变量多而目标少时(形状优化通常如此)伴随路线便宜得多;设计变量少而目标多时,前向路线更合适。
This page has and : the tangent route solves 9 linear systems of order 9, 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.
设计导数与反问题
Design derivatives and inverse design
一次伴随求解之后,每个导数只是内积
After one adjoint solve, each derivative is an inner product
前向模式要把设计种子依次取为单位向量 ,才得到梯度的各个分量;反向模式一次给出全部分量 。本例源项导数非常稀疏:格心一个源区间影响一个控制体;格点一个源区间影响左右两个半控制体,但入口约束行不保留源项。伴随把这些局部直接影响与全局目标权重配对。
Forward mode needs the design seed for each in turn to obtain the gradient components; reverse mode gives all of them, , at once. Source derivatives are sparse: one cell interval affects one volume; one node-based interval affects two half-volumes, with the inlet constraint row excluding its source contribution. The adjoint pairs these local direct effects with global objective weights.
负梯度给出目标局部下降方向,但大步长可能越过线性近似有效的范围。演示中的 BFGS 利用历史梯度差构造逆 Hessian 近似,再用回溯缩短步长,直到目标充分下降。每个候选设计都需要重新解原始方程;伴随提供的是该设计处的梯度。
The negative gradient is a local descent direction, but a large step can leave the range where linear prediction is useful. The demo’s BFGS uses gradient differences to approximate an inverse Hessian, then backtracks until the objective decreases sufficiently. Each candidate design requires a fresh primal solution; its adjoint supplies the gradient there.
一次伴随求解即可计算所有源项和入口导数;本例目标函数对设计变量没有显式依赖,所以梯度只剩伴随与残差设计偏导的内积,也就是反向模式的设计伴随 。tangent 对每个设计种子 单独求解,结果应与这个内积一致。
One adjoint solve yields the derivatives for every source and for the inlet. The objective has no explicit design dependence here, so the gradient is just the adjoint contracted with the residual design derivatives: the design adjoint of reverse mode. A tangent solve handles one design seed at a time and must agree with that contraction.
源项导数必须清零入口约束行;入口值的直接残差导数只有 ,所以 。源项首区间仅对剩余半控制体贡献 ,这会影响入口局部灵敏度的一致性。
Source derivatives must zero the inlet constraint row. The inlet’s direct residual derivative is only , so . The first source interval contributes only to the retained half-volume, affecting local inlet-sensitivity consistency.
设计伴随从哪里来
Where the design adjoint comes from
以 为种子运行第 9 节的反向代码,并让它同时累加设计输入的伴随。读取设计的只有两类语句,所以 只从两处收到贡献。源项来自反向语句 7:,再按 对 的偏导分给源区间,;第 0 行的种子已被语句 8 的伴随清零,所以 只从节点 1 得到贡献。入口值来自反向语句 8:。合起来:
Run the reverse code of Section 9 with the seed and let it also accumulate the adjoints of the design inputs. The design is read by two kinds of statement only, so collects from two places. The sources come from reverse statement 7: , passed on to the source intervals by the partial derivatives of with respect to , ; the seed of row 0 has been cleared by the adjoint of statement 8, so receives only the contribution of node 1; the inlet value comes from reverse statement 8: . Together:
图 16把这一次反向扫描逐句播放:种子 从残差端出发,状态伴随 与设计伴随 用 += 逐步累加,最后 (伴随方程成立), 就是梯度。图 17把 与梯度的每个分量画成一张灵敏度图,并用中心差分逐个核对;与图 14的伴随结果一样,这些分量只用了一次伴随求解。
Figure 16 plays this reverse sweep statement by statement: the seed starts at the residual end, the state adjoint and the design adjoint accumulate with +=, and at the end (the adjoint equation holds) and is the gradient. Figure 17 draws and every gradient component as a sensitivity map and checks each by a central difference; like the adjoint results of Figure 14, these components cost a single adjoint solve.
用梯度调整源项
Use the gradient to adjust the source
反问题演示只优化源项,固定入口,避免入口与首源项近似互相补偿。每次 BFGS 步都重新收敛原始方程并计算伴随,使用 Armijo 回溯接受下降步;重置或切换格式会恢复基准设计。
The inverse-design demo optimises only sources and fixes the inlet, avoiding near compensation between the inlet and first source. Each BFGS step reconverges the primal and adjoint equations and uses Armijo backtracking to accept a descent step. Resetting or changing discretisation restores the base design.
有限差分:检查总导数
Finite differences: check total derivatives
用不经过伴随的计算作对照
Check with a computation that does not use the adjoint
有限差分直接改变设计、重新求解、再比较目标。步长 是一个数(不是网格宽度 ,也不是求导记号)。对 和 分别作 Taylor 展开,相减后偶次项消去;除以 得到中心差商 ,它的首个截断误差是二阶。
Finite differences change the design, solve again, then compare objectives. The step is a number (neither the mesh width nor a derivative symbol). Taylor-expand and ; subtraction cancels the even powers, and division by gives the central quotient , whose leading truncation error is second order.
但两次目标都有舍入和求解误差,作差再除以很小的 会把这些误差放大。因此应看一组步长:先出现二阶下降区,再出现平台或回升。如果整段都不接近伴随,先检查原始残差是否收敛,以及两边是否使用同一个固定目标。
Both objective values also contain rounding and solve error; subtracting and dividing by tiny amplifies it. Inspect a range of steps: a second-order region, then a plateau or rise. If none approaches the adjoint result, first check primal convergence and that both solves use the same frozen target.
中心差分的完整求解流程
Complete central-difference workflow
freeze target_field
for tau in decreasing_step_sizes:
u_plus = solve_primal(D + tau*e[k])
u_minus = solve_primal(D - tau*e[k])
FD_c = (J(u_plus) - J(u_minus)) / (2*tau)
compare FD_c with adjoint_gradient[k]固定目标场,对当前设计做正负扰动,每次重新收敛完整原始残差。中心差分误差先按 减小,随后受舍入误差与求解容差影响。不能仅根据最小步长的结果判断伴随是否正确。
Freeze the target field, perturb the current design in both directions, and reconverge the full primal residual each time. Central-difference error initially decreases as , then is affected by roundoff and solve tolerance. The smallest step alone is not a reliable adjoint check.
前向差分截断误差通常为 ,中心差分通常为 ;两者都会在小步长下放大目标值的舍入和求解误差。这里的步长 是一个数,即设计变量的差分步长;它不是网格宽度 ,也不是求导记号。目标场始终固定,不随着扰动重新生成。
Forward differences usually have truncation error, central differences . Both amplify objective roundoff and solve errors at small steps. Here the step is a number, the finite-difference step of the design; it is neither the mesh width nor a derivative symbol. The target remains frozen throughout the check.
下面可独立选择任一源区间或入口值,也可取全部变量的最大误差。扫描覆盖 到 ,可比较前向和中心差分,并输入其他步长查看各分量的目标值及收敛残差。切换空间精度、网格或边界后,所有检查针对新的完整离散系统重新计算;切换语言或差分显示复用当前模型的缓存。
Choose any source interval or the inlet value locally, or inspect the maximum error over all variables. The sweep covers through ; compare forward and central differences and enter another step to inspect objectives and convergence residuals. Changing spatial order, mesh or boundary rebuilds checks for the new full discrete system. Language and display changes reuse the current model cache.
读图时先选入口、中部源项和末端源项,比较误差最低点是否相同,再切换全部变量。中心差分不一定在每个步长都更准确;误差已经进入平台时继续减小步长,可能使误差反而增加。本面板缓存基准、正扰动、负扰动三组解,方便直接切换两种差分公式。
Compare the inlet, a middle source and the final source before selecting all variables: their minimum-error steps can differ. Central differences need not be more accurate at every step. Reducing an already tiny step can increase error. This panel caches baseline, positive and negative perturbed states so either formula is immediately available.
四种误差要分开:空间离散误差由网格加密检查;非线性求解误差由原始残差检查;线性求解误差由 tangent/adjoint 真残差检查;求导实现由转置、差分与复数检查。点积相等只证明两条传播互为转置;共同遗漏同一项也可能通过。近零梯度分量请同时查看绝对误差。
Separate four errors: spatial error through refinement; nonlinear solve error through primal residuals; linear solve error through true tangent/adjoint residuals; derivative implementation through transpose, finite differences and complex checks. A dot test alone only establishes paired transposes; jointly omitting a term may still pass. Inspect absolute errors for near-zero gradient components.
Taylor 实验使用固定设计方向(种子),对每个步长 重新收敛原始问题,并保持目标场不变。,,其中 是伴随梯度给出的方向导数。光滑且 时,预期 为一阶、 为二阶;驻点或特殊方向可能有更高阶。步长过小后,求解误差与浮点消减会破坏这些斜率,不能把这一段当作空间精度。
The Taylor experiment reconverges the primal for each step along a fixed design direction (seed) with a fixed target. and , where is the directional derivative from the adjoint gradient. For a smooth problem and nonzero , expect first- and second-order behaviour; stationary points or special directions can give higher orders. At tiny steps, solve errors and cancellation spoil the slopes. These slopes are not spatial accuracy.
复数步长:用虚部检验导数
Complex step: check derivatives with an imaginary perturbation
为什么一个很小的虚部能携带导数
How a tiny imaginary part carries a derivative
把普通 Taylor 展开中的实步长换成 。对于可解析延拓的函数,虚部首项恰好是 ,再除以 就得到导数。这个过程没有把两个接近的实数相减,因此能使用很小的步长。
Replace the real Taylor step by . For an analytic extension, the leading imaginary part is ; divide by to obtain the derivative. There is no subtraction of nearly equal real values, allowing a very small step.
复数目标必须延续原来的代数平方 ,不能改成带共轭的模平方 ;后者不是这里需要的解析延拓,会抹去虚部导数。绝对值、比较分支和限制器也需要单独分析。
The complex objective must retain the algebraic square , not the conjugate modulus square . The latter is not the required analytic extension and removes the imaginary derivative. Absolute values, comparisons and limiters require separate analysis.
重构、固定耗散通量和目标函数在本例中均为解析运算。将设计扰动为 (步长 是一个实数),保持目标场为实数常量,用独立复数残差收敛复数状态。虚部残差按 缩放后也必须收敛;只检查实部会错误地提前停止。
Reconstruction, the fixed-dissipation flux and the objective are analytic here. Perturb the design to (the step is a real number), keep the target real and fixed, and converge a complex state using an independent complex residual. The imaginary residual divided by must also converge; checking only the real part would stop too early.
复数步长验证使用冻结在收敛实解处的实数 GS 预处理,不经过 GMRES 的范数和 Arnoldi 分支。它验证收敛离散方程的导数;不声称 GMRES 与复数迭代逐步相同。本例没有 minmod、绝对值通量开关或限幅分支;加入它们后必须重新分析可微性。
The complex-step check uses a real GS preconditioner frozen at the converged real state, bypassing GMRES norms and Arnoldi branches. It validates derivatives of the converged discrete equations; it does not claim iteration-by-iteration equality with GMRES. There are no minmod, absolute-value flux switches or clipping branches here; adding them requires a new differentiability analysis.
图 18把第 13 节的差分步长扫描延伸到复数步长。取源项 (第 13、14 节面板的默认分量),格式同页首默认(,一阶)。所有目标值都来自本页的重新求解:前向与中心差分用第 13 节同一个原始求解器,每个扰动都重新收敛;复数步长用本节的复数迭代,固定做 60 次,远超停止准则,让缩放残差停在舍入平台上。误差都相对于直接法求出的伴随梯度。
Figure 18 extends the step sweep of Section 13 to the complex step. It uses the source (the default component of the Section 13 and 14 panels) and the default discretisation (, first order). Every objective value comes from this page's re-solves: the forward and central differences use the primal solver of Section 13, reconverged for every perturbation; the complex step uses this section's complex iteration, run for a fixed 60 iterations, far past the stopping criterion, so that the scaled residual rests on its round-off floor. Errors are relative to the adjoint gradient from a direct solve.
前向和中心差分的误差先按斜率 1 和 2 下降,分别在 和 附近到达 V 形底部,然后回升:步长太小时,两个几乎相等的 相减,原始求解的收敛容差和舍入误差被 放大。复数步长没有这一步相减。它在大步长时与中心差分重合,因为两者的首项截断误差都是 (符号相反); 以后停在 左右的平台,与 同一量级,一直保持到 。
The forward and central errors first fall with slopes 1 and 2, reach the bottom of the V near and , and then rise: at tiny steps two nearly equal values of are subtracted, and the convergence tolerance of the primal solve and round-off are amplified by . The complex step makes no such subtraction. At large it coincides with the central difference, because both have the leading truncation error (with opposite signs); for it settles on a plateau near , the order of , and stays there down to .
点积检验通过了,为什么还要检查残差
Why a dot test is not enough
故障对照分别检查转置、残差微分与总导数。二阶时把两个方向都替换为一阶 Jacobian,点积检验仍然成立,但复数残差导数检验会失败;这是“自洽但求错导数”的例子。
The fault comparison checks transposition, residual differentiation and total derivatives separately. At second order, replacing both operators by a first-order Jacobian still passes a dot test, but fails the complex-residual derivative check: the pair is internally consistent yet differentiates the wrong residual.
复数步长求解:实部和缩放虚部同时收敛
Complex-step solve: converge real and scaled imaginary residuals
u_real = solve_primal(D); u_imag = zeros(M)
tau = 1e-30; D_imag = tau*e[k]
B = fixed_GS_preconditioner_at(u_real)
repeat:
(r_real, r_imag) = complex_residual(u_real, u_imag, D, D_imag)
if max(norm(r_real), norm(r_imag)/tau) <= tolerance: break
u_real -= B(r_real)
u_imag -= B(r_imag)
CS = Σ W[i]*(u_real[i] - target[i])*u_imag[i]/tau
# Target stays real and fixed; the objective uses an algebraic square.网格收敛:一阶与二阶
Mesh convergence: first versus second order
本节实验:线性重构精确、PDE 解达到二阶、目标或梯度达到二阶,是三种不同的结论。分别检查重构、加密解和方向梯度。
Experiment: Exact linear reconstruction, second-order PDE solutions and second-order functionals or gradients are distinct claims. Check reconstruction, refined solutions and directional gradients separately.
先收敛迭代,再比较离散误差
Converge the iteration before measuring discretisation error
同一网格上把残差降到很小,只说明代数方程被解准了;它不说明这个离散解离连续解有多近。空间精度需要在越来越细的网格上重复求解。若 ,把 减半,误差应约缩小为原来的 :一阶约一半,二阶约四分之一。
A small residual on one mesh means the algebraic equations are solved accurately. It does not measure distance to the continuous solution. Spatial accuracy requires repeated solves on finer meshes. If , halving reduces error by about : one-half for first order, one-quarter for second order.
格心值是平均值,所以参考值 也必须是控制体上的精确平均,而不是直接拿中心点精确值来比。两者对光滑函数本身相差 ,在研究二阶时不能忽略。格点值则与对应点值比较。
A cell unknown is an average, so its reference must be the exact volume average, not simply the exact centre value. These differ by for smooth functions and cannot be interchanged in a second-order study. Node values are compared with point values.
把区间逐次二等分,在 上分别求解,就能检查一阶与二阶格式的实际精度。格心未知量与精确单元平均值比较,格点未知量与精确节点值比较。采用控制体加权 误差,并报告相邻网格观测阶。
Halve the grid spacing successively and solve on to measure the accuracy of both discretisations. Cell unknowns are compared with exact cell averages; node unknowns with exact point values. We report control-volume-weighted errors and observed orders between successive grids.
tangent 加密使用固定的物理设计方向:每个源区间的种子 (源项密度处处增加 1),入口种子 ;对应的连续切向量是 (附录 A 用连续函数的变分写同一个方向)。不能在不同网格上只给一个宽度随网格变化的源区间加种子,却把它称为同一个连续方向。目标误差使用共同的连续目标场投影和连续积分 ;它与第 7 节用于反问题的离散目标定义有意区分。
Tangent refinement uses a fixed physical design direction: the seed on every source interval (the source density raised by one everywhere) and the inlet seed . The corresponding continuous tangent is (Appendix A writes the same direction as a variation of continuous functions). Seeding a single source interval whose width changes with the mesh would not represent the same continuous direction. Objective error uses a shared projected continuous target and a continuous integral for , deliberately distinct from the discrete inverse-design target in Section 7.
空间误差图的横轴是网格宽度;求解器收敛图的横轴是迭代或工作量。这两种“收敛”不能混用。伴随全域误差还受边界闭合影响,第 16 节单独讨论。
Spatial-error plots use mesh width on the horizontal axis, whereas solver convergence uses iteration or work. These are different notions of convergence. Boundary closure also affects global adjoint errors, treated separately in Section 16.
伴随一致性:结论与阅读范围
Adjoint consistency: conclusions and scope
连续终端值问题在当前正值光滑解上适定。离散伴随的内部收敛与全域收敛可以不同:入口闭合造成的局部边界层,不能当作连续问题不适定的证据。格点入口约束乘子须与场值分开比较。
The continuous terminal-value problem is well posed on the current smooth positive branch. Interior and global discrete-adjoint convergence can differ: a local layer caused by inlet closure does not prove continuous ill-posedness. Compare the node inlet constraint multiplier separately from field samples.
点积测试检查离散转置,差分和复数步长检查同一离散目标的梯度。独立连续参考回答另一问题:加密后是否趋向相同的连续灵敏度。误差不下降时,先核对边界、目标、乘子含义、正则性和参考误差,再讨论适定性。
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.
图 19在整个区间上比较离散伴随与连续伴随(一阶,)。远离入口, 随加密逼近 : 处的最大误差依次为 、、,最细两套网格之间的观测阶为 ,接近一阶格式的 1。入口第一个控制体的误差却不趋于零,而是趋向 (依次为 、、)。这是已知的入口伴随边界层:入口处的离散伴随方程缺少上游的一项,相当于入口外的伴随值取零。它不影响光滑方向导数的收敛,却使全域 误差只有半阶。附录 A 推导了这个极限,见图 A.3与图 A.5。
Figure 19 compares the discrete and the continuous adjoint over the whole domain (first order, ). Away from the inlet, approaches under refinement: the largest error over is , , , with observed order between the two finer meshes, close to the 1 of a first-order scheme. The error in the first control volume at the inlet does not go to zero; it tends to (here , , ). This is the known inlet adjoint boundary layer: the discrete adjoint equation at the inlet lacks the upstream term, as if the adjoint outside the inlet were zero. It does not spoil the convergence of derivatives in smooth directions, but it limits the global error to half order. Appendix A derives the limit; see Figure A.3 and Figure A.5.
用伴随估计目标误差
Estimate objective error with the adjoint
伴随还可把细网格残差转换为目标误差估计。将粗解线性插值到细网格( 个区间)得到 ,在线性化点求解 ,则 。这个估计针对同一细网格目标,不把不同网格的离散目标直接相减当作误差。
The adjoint also converts a fine-grid residual into an objective-error estimate. Interpolate the coarse solution linearly onto the fine grid ( intervals) to obtain , and solve . Then . This estimate concerns the same fine-grid objective; it does not directly subtract objectives with different discrete targets.
离散乘子与连续场为什么不能直接混用
Why a discrete multiplier is not automatically a continuous field
正文使用普通数组点积定义 。若另选带权内积,算子的伴随就会改变:状态权矩阵为 、残差权矩阵为 时,带权伴随是 。代码的普通转置不能在没有换元说明的情况下被写成该带权算子。
The main derivation defines using ordinary array dot products. A weighted inner product changes the adjoint: with state and residual weights and it is . The coded ordinary transpose cannot be relabelled as that weighted operator without a change of variables.
残差的归一化也会改变乘子的数值。若改为 , 为可逆的行尺度矩阵(例如各控制体的体积),则在收敛解上 ,;目标总梯度保持不变。即使 随设计改变,由 产生的项也乘以 ,在精确收敛时为零。把控制体积分残差改成体积平均残差后,不能仍直接比较原来的乘子曲线。
Residual normalization also changes multiplier values. For with an invertible row-scaling matrix (for example the control-volume sizes), at convergence and , while the total objective gradient is unchanged. Even when depends on the design, the terms from multiply and vanish at exact convergence. Multipliers of integrated and volume-normalized residuals cannot be compared without rescaling.
连续附录从微分方程、积分目标和分部积分出发;正文从实际离散残差及其有限维链式法则出发。两者可以比较同一设计变量的目标总导数与加密极限,但还须匹配符号、变量、积分权重和边界约束。被替换的边界行的乘子尤其不能直接解释为连续场采样值。
The continuous appendix starts from differential equations, integral objectives and integration by parts; the main text starts from actual discrete residuals and finite-dimensional chain rules. Compare total derivatives of the same objective and design variable and their refinement limits, after matching signs, variables, weights and boundary constraints. Multipliers of replaced boundary rows especially are not direct samples of a continuous field.
GS 固定点与 GS 预处理 GMRES
GS fixed point and GS-preconditioned GMRES
从固定方向修正到组合多个方向
From repeated corrections to combining search directions
把前向或伴随系统统一写成 。GS 固定点每次采用 作为修正。若误差中有衰减很慢的分量,就会需要很多轮。GMRES 保存已经产生的方向,在它们张成的子空间中选择使残差二范数最小的组合;GS 则先把每个新方向变成较容易被 处理的预处理方向。
Write either linear system as . GS fixed point repeatedly corrects by . Slowly decaying error components can require many iterations. GMRES retains directions and chooses a combination minimising the residual norm in their span. GS first maps each new direction to a preconditioned direction for the action of .
Arnoldi 正交化逐步建立正交基 ,把大矩阵在这组方向上的作用压缩成小的 Hessenberg 矩阵。对小矩阵解最小二乘问题,就能选出组合系数。右预处理时实际更新使用 ;保留真实残差 ,便能与固定点方法按同一标准比较。
Arnoldi builds an orthonormal basis and represents the large operator on those directions by a small Hessenberg matrix. A small least-squares solve supplies the combination coefficients. With right preconditioning, the state update uses . Monitoring the true residual allows comparison with fixed point on the same basis.
固定点方法反复沿 GS 给出的修正方向更新。GMRES 则在 Krylov 子空间中组合方向,求使残差二范数最小的近似解。GMRES 仍是迭代法,不是直接分解求解器,也不保证每个小问题都比 GS 更省工作。
Fixed point repeatedly updates in the correction direction supplied by GS. GMRES combines directions in a Krylov subspace to minimise the residual 2-norm. It is still an iterative method, not a direct factorisation solver, and need not use less work than GS on every small problem.
让 GS 做预处理
Use GS as the preconditioner
采用右预处理, 是固定次数、零初值的 GS 扫描;精确算子始终是 或 ,一阶 Jacobian 仅用于 GS 预处理。伴随使用对应反向扫描 。下式的 是 Arnoldi 过程得到的 上 Hessenberg 矩阵。
Right preconditioning uses , a fixed number of zero-initialised GS sweeps. The exact operator remains or ; the first-order Jacobian is used only in the GS preconditioner. The adjoint uses the corresponding reverse sweeps . Below, is the upper Hessenberg matrix built by the Arnoldi process.
每一步先计算 ,再计算 ,用两遍修正 Gram–Schmidt 正交化,Givens 旋转更新小型最小二乘问题。达到重启长度后,用当前真实残差重新建立子空间。页面记录的是 ,不是预处理残差。
Each step applies , then , orthogonalises with two-pass modified Gram–Schmidt, and updates the small least-squares problem with Givens rotations. At the restart length, rebuild the subspace from the current true residual. The page records , not the preconditioned residual.
每次真实残差检查也计入矩阵向量乘次数。比较同时列出 GS 扫描次数与 Krylov 存储量。本实现固定预处理;若在一次线性求解内改变 GS 扫描次数或用变化的内层容差,应使用适合可变预处理的 FGMRES。
True-residual checks are included in matrix-vector-product counts. Comparisons also report GS sweep counts and Krylov storage. This implementation has fixed preconditioning; changing GS sweep counts or using a variable inner stopping tolerance during one linear solve calls for a method such as FGMRES.
重启右预处理 GMRES:与本页计算对应
Restarted right-preconditioned GMRES used on this page
K = exact_forward_or_adjoint_product
B = fixed_zero_start_GS_or_reverse_GS
x = zeros(M)
while norm(b - K(x)) > tolerance:
x_base = x; r = b - K(x); beta = norm(r)
V[0] = r / beta; g_small = beta*e[0]
for j = 0,...,min(M,restart_length)-1:
Z[j] = B(V[j]); w = K(Z[j])
orthogonalise w against V[0:j] twice; store H[:,j]
h_next = norm(w)
if h_next > breakdown_threshold: V[j+1] = w/h_next
apply previous Givens rotations to H[:,j]
rotate H[j:j+1,j] and g_small[j:j+1] with a new rotation
y = back_substitute(H[0:j,0:j], g_small[0:j])
x = x_base + Σ Z[i]*y[i]
if norm(b - K(x)) <= tolerance: return x
if h_next <= breakdown_threshold: break
# Restart with the true residual; enforce the total iteration limit.比较工作量时,不能只比较迭代编号:GMRES 还做正交化、内积并保存基向量。页面同时列出精确矩阵乘积次数、GS 扫描次数、存储量与真实残差;重启长度是内存和收敛速度之间的一个选择。
Iteration counts alone do not compare work: GMRES also orthogonalises, takes inner products and stores basis vectors. The page reports exact products, GS sweeps, storage and true residuals. Restart length trades memory against convergence.
总结:同一套网格、通量与重构,三种传递方向
Summary: one mesh, flux and reconstruction, three directions of data flow
把前面的计算并排放在一起,可以看到:原始、前向和伴随使用同一张网格、同一个通量和同一套重构系数。改变的是传递的数据,以及收集与分发的方向。
Putting the calculations side by side shows that the primal, forward and adjoint share the mesh, flux and reconstruction coefficients. What changes is the data being propagated and the direction of gather and scatter.
图 21把这三种计算并排画出:方框是存储的数组,箭头是数据流,颜色表示角色。
Figure 21 draws the three calculations side by side: boxes are stored arrays, arrows are data flow, colours mark the role.
| 计算Calculation | 原始Primal | 前向(tangent)Forward (tangent) | 伴随(adjoint)Adjoint |
|---|---|---|---|
| 结果Result | |||
| 面上数据Data on a face | |||
| 重构的作用Role of reconstruction | 由未知量重构面值Unknowns to reconstructed face states | 用相同系数组合切向量Combine tangents with the same weights | 按转置权重累加(+=)到模板各点Accumulate (+=) to stencil points with transposed weights |
| 要求解的方程Equation to solve | |||
| GS 修正GS correction |
- 空间精度。一阶在控制体内取常值;二阶把斜率带到面上。前向与伴随都要经过这一步重构,边界模板也包括在内。
- Spatial accuracy. First order uses constant states; second order carries slopes to the faces. Both derivatives include this reconstruction and its boundary stencils.
- 求解。一阶 Jacobian 与 GS 扫描给出修正。固定点迭代反复使用这个修正;GMRES 把它作为预处理,在 Krylov 子空间中组合方向。两者都以完整方程的残差判断收敛。
- Solving. The first-order Jacobian and GS sweeps supply corrections. Fixed point repeatedly applies them; GMRES uses them as a preconditioner and combines directions in a Krylov subspace. Both check the residual of the full equation.
- 梯度。前向求解给出一个设计种子 对应的状态切向量 ;一次伴随求解可以计算所有源项和入口值对目标函数的影响,即设计伴随 。
- Gradients. A forward solve gives the state tangent for one design seed . One adjoint solve gives objective sensitivities to every source and to the inlet value, the design adjoint .
- 验证。点积检查转置关系,有限差分与复数步长检查导数,网格加密检查空间精度。它们各自回答不同的问题。
- Verification. Dot products check transposition, finite differences and complex step check derivatives, and mesh refinement checks spatial accuracy. Each addresses a different question.
- 连续伴随。本例正解分支上的连续伴随适定;离散伴随能否逼近它,还取决于边界处理及乘子的含义。内部误差、全域误差与入口灵敏度需要分别比较。
- Continuous adjoint. The continuous adjoint is well posed on the positive branch. Its discrete approximation also depends on boundary treatment and the interpretation of multipliers. Interior error, global error and inlet sensitivity are compared separately.
交互面板中的解、残差、前向与伴随、梯度都由页面实时计算。网格加密图使用同一模型预先求出的数据。一阶/二阶、边界和语言切换会更新相应说明与结果。
Solutions, residuals, forward and adjoint products, and gradients in the interactive panels are computed live. Refinement plots use results computed in advance with the same model. Spatial-order, boundary and language controls update the corresponding explanations and results.
附录 A:连续伴随在简单算例中的验证与对比
Appendix A: verifying and comparing continuous adjoints on simple cases
本项目以离散伴随为主线:从可逐句检查的小模型,逐步走向工业级非结构网格 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.
连续 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
A.1 完整推导与适定性
A.1 Full derivation and well-posedness
以下第 1、2 部分摘自全项目通用的连续伴随说明,其中提到的二维物体、对流扩散和拟一维 Euler 算例见对应页面;第 3 部分是本页的一维 Burgers 推导。
1. 符号、积分测度与三类不同的检查
将定常 PDE 写为 ,目标泛函记为 ,采用拉格朗日泛函
页面求解的是 。这里的连续乘子采用相反符号:在守恒方程上有 。这些残差已经是控制体积分后的守恒平衡,不需要再除以控制体的体积或面积。被代数边界约束替换的行对应约束乘子,其数值不是连续场的点值。
三类检查分别回答不同的问题:
- 转置/点积测试检查 的实现。
- 对同一离散问题做有限差分,检查 。
- 将独立求得的连续伴随与一系列加密网格上的离散伴随比较,检查指定范数下的一致性。比较时必须匹配 PDE、边界条件、目标泛函、计算区域与目标数据。
有限组网格上误差下降是数值证据,不是收敛定理。即使积分范数下的误差下降,边界层或特征跳跃仍可能阻止最大范数下的收敛。
2. 适定性与伴随一致性不是同一件事
出现差异提供了诊断机会,但差异本身不能证明问题不适定。应先明确状态空间与边界迹空间、允许的数据、目标泛函以及范数。Hadamard 意义下的适定性要求解存在、唯一且连续依赖于数据。正确推导的连续伴随可以是适定的,而离散边界闭合却不一致;一个边界数据与目标泛函不匹配的适定方程,仍然是错误的伴随。约束乘子、有限正则性和参考解误差也会造成表观差异。
线性化、零空间与稳定性
把原始边界条件的齐次线性化约束纳入算子 的定义域 。如果该算子是适当空间之间的有界同构,则其 Banach 伴随也是相应对偶空间之间的有界同构。这给出一种充分的局部分析框架,并不是写出形式微分表达式就能得到的结论。对于无界 PDE 算子,必须先明确其定义域及图范数,或者采用适当的弱形式。对于非线性原始问题,仅有存在唯一性还不能保证解映射可微;需要在所选解支上检查局部光滑性及线性化可逆性。
对于 ,原始算子的零模态给出必要的兼容性条件:
若 含有非零元素,则伴随一旦有解就不唯一。仅有与原始零空间正交的条件,并不能构成一般的无穷维存在性定理:还需要考虑值域闭性和稳定性。对齐次伴随边界数据,所需估计具有以下形式:
如果边界数据非齐次,还应加上相应的边界迹范数。每个网格尺度下逆算子有界也不够:在正确加权的离散范数下,稳定界必须与网格尺度无关。不考虑积分残差和控制体权重的原始欧氏条件数,本身不能作为连续问题的诊断结论。
一个可显式求解的对流例子
取 ,固定原始入口值,目标变分为 。对原始线性化 使用 Green 恒等式,得到
- 缺少出口数据: 解可以加上任意常数,因此不唯一。
- 额外照搬原始入口条件 : 存在性要求 。当 时,正确解是 ;额外入口条件使问题无解。
- 失去一致稳定性: 有 。取 就达到 的放大倍数。每个固定的 都可以对应适定问题,但对于这类数据,不能得到在 时仍统一有效的界。速度为零时,原始算子退化。
在原二维物体算例中,法向指向流体计算域的外部,所以迎风物面是流体域的原始出流边界。一条已经携带 的特征线,不可能在这里同时满足零物理通量。首先失败的是所假设的光滑正值原始问题,而不只是某个已经建立的连续伴随的离散逼近。在假定的壁面状态 处,
对任意变分 都成立:通量约束丢失了一阶约束信息,特征速度也同时退化。把壁面条件照抄到形式伴随上,不能恢复有效的光滑解支。弱解/熵解表述需要另行分析边界迹与可微性。
对于特征边界对照,每条内部射线都会在有限路程内到达出流边界。边界数据几乎处处确定内部场;切线或角点引出的跳跃可以在允许非光滑解的空间中处理,本身不意味着不适定。特征边界附近的比值 需要兼容性分析:在这里的物面目标中,分子与分母具有可相消的法向因子,不能只看分母就断言伴随发散。
其余算例说明了什么
对流扩散。 正扩散系数提供椭圆主部,但这不是任意低阶项与 Robin 边界系数下的普遍可逆性证明。仍需检查齐次零空间,以及所用能量估计或最大值原理的假设。独立计算支持当前参数组,没有证明对所有参数都成立的定理。把目标所要求的壁面伴随值 1 改成 0,可能得到另一个适定问题,但会给出错误的目标灵敏度。在远场同时规定状态值和法向导数,一般会过约束。把被替换方程行的乘子当成壁面场值,则是另一类代数含义错误。
拟一维 Euler。 对所述正向亚声速分支,特征速度 解释了原始入口两个条件、出口一个条件,以及伴随的相反条件数。条件数目正确是必要条件,不是充分条件:特征传播与边界约束共同组成的边值算子还必须可逆。声速点处通量 Jacobian 奇异,当前使用 的表示失效,需要喉部兼容性/正则性分析。这不意味着所有跨声速伴随都不存在。激波还需要考虑激波位置变分及相应内部条件。这些现象超出当前数值证据的范围,可参考 Giles 与 Pierce 的连续喷管伴随分析。强、弱边界处理对数值伴随的不同影响,见 Duivesteijn 等人的研究。
因此,教学顺序应当是:先明确连续原始问题及所选解支;推导伴随,包括其边界定义域;分析存在性、唯一性和稳定性;随后才检查对该问题的一致、稳定逼近。反例应明确保留其失败的前提,不能悄悄替换模型后宣布验证成功。
3. 一维 Burgers:二阶重构、连续伴随与边界一致性
光滑正解分支与连续 tangent
考虑固定区间 [0,1] 上的定常方程,入口值为正。令源项的原函数从入口积分,选取正平方根:
只要根号内严格正且远离零,光滑解唯一,源项和入口扰动对解的映射稳定。令 v 为状态扰动,线性化得到:
教学基准采用 s=0.5、入口值 1,故 u(x)=sqrt(1+x)。目标场对应 s=sin(πx),在求导期间固定。网格加密的 tangent 参考使用同一个物理扰动 δs=1、=0.2,而不是随网格缩小的单个源项区间。
从变分到终端伴随方程
目标函数是固定目标场的平方误差积分。令连续伴随 λ 满足以下终端值问题:
这里是 -uλ′,不是 -(uλ)′。原因是 tangent 微分算子作用于乘积 uv;其形式转置由下面的积分分部直接确定:
代入 tangent 方程即可消去 v,得到连续总导数及可直接计算的伴随积分:
于是区间常源项的连续灵敏度是该区间内 λ 的积分,入口灵敏度是 λ(0)。页面使用复合五点 Gauss 积分计算参考值,独立于离散伴随线性求解。
适定性与离散不一致是不同问题
若 u≥c>0,则终端条件唯一确定 λ,且有稳定性估计:
因此当前基准的连续伴随适定。离散伴随出现入口边界层,不应被解释为这个连续终端值问题不适定。若 u 接近或穿过零,以上稳定性界不再统一成立;若有激波,还需要弱解、激波位置变化及相应的伴随传递条件,不能套用这里的光滑推导。
二阶重构不会自动带来全域二阶伴随
格心未知量是单元平均值,格点未知量是节点值。二阶采用内部中心斜率、端点单边斜率,重构到真实面位置后计算通量。格心出口必须先重构到 x=1;格点第一内部对偶面位于 h/2,迎风通量也必须使用该位置的重构迹。
用 、 表示重构,用 B 表示有向面分发,忽略此行公式中单列的边界常数及约束替换,则精确二阶导数为:
页面使用离散拉格朗日函数 J+ψTR,因此 A₂Tψ=-,平衡方程对应的连续场应比较 λ 与 -ψ。格点第 0 行是强约束,乘子 ψ₀ 不能直接当作控制体内 λ 的取值;改为比较入口灵敏度 -ψ₀ 与 λ(0)。
入口约束会移除该行的源项贡献;重构的转置又会把其他残差行的权重传回边界附近未知量。因此必须同时检查内部场、全域场、入口导数和首源项导数。当前计算显示:默认耗散入口下,二阶内部误差接近二阶,但全域误差受有限幅值局部边界层影响,接近半阶。格心迎风入口改善全域收敛;格点二阶迎风对照仍存在边界层。
首源项区间随网格缩小,其导数本身会趋于零;绝对差趋零不足以证明局部一致性。因此页面额外报告离散导数与连续区间积分的比值。对固定光滑物理扰动,则另外报告总方向导数误差。
GS 与 GMRES 只改变求解过程
令 P=A₁(u*) 是在收敛二阶状态处计算的一阶近似 Jacobian。每次从零开始执行固定次数 GS 扫描,形成线性近似逆 。tangent 与 adjoint 的固定点迭代使用精确残差,GMRES 也使用精确二阶算子:
一次前向 GS 对应 B₁=(diag(P)+L)⁻¹;伴随采用对应的反向扫描 B₁T=(diag(P)+LT)⁻¹。这里的反向 GS 解的是 PT 系统。多次扫描也固定次数、零初值,以保持标准 GMRES 的预处理线性且固定。换成 GMRES 或把线性残差收敛得更小,都不能消除由离散边界闭合造成的伴随不一致。
Parts 1 and 2 below come from the project-wide notes on continuous adjoints; the 2-D body, advection–diffusion and quasi-1-D Euler cases they mention have their own pages. Part 3 is the one-dimensional Burgers derivation for this page.
1. Signs, measures and three distinct checks
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 transpose/dot test checks the implementation of .
- Finite differences of the same discrete problem check .
- Comparing an independently solved continuous dual with a sequence of refined discrete duals tests consistency in a stated norm. It requires matching PDE, boundary conditions, functional, domain, and target data.
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.
2. Well-posedness is not the same as adjoint consistency
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.
Linearization, kernels and stability
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.
An explicit transport example
For , fixed primal inflow and objective variation , the Green identity for the primal linearization gives
- Missing outlet data: an arbitrary additive constant remains, so uniqueness fails.
- An extra copied primal inlet datum : existence requires . With , the correct solution is ; the extra inlet condition makes the problem unsolvable.
- Loss of uniform stability: . Taking attains the amplification. Each fixed is well posed, but there is no bound uniform as for this data class. At zero speed the original operator degenerates.
In the original 2-D body example, normals point out of the fluid domain: the upwind body is therefore primal outflow. A characteristic already carrying cannot satisfy zero physical flux there. This is first a failure of the assumed smooth positive primal problem, not merely a poor discrete approximation to an otherwise established dual. At a putative wall state ,
for every variation : the flux constraint loses its first-order information while the characteristic speed degenerates. Copying the wall condition onto a formal dual cannot restore a valid smooth solution branch. Weak/entropy formulations would need separate trace and differentiability analysis.
For the characteristic-boundary contrast, every interior ray reaches outflow in finite distance. Boundary data determine the field almost everywhere; tangency/corner jumps are compatible with a nonsmooth solution space and do not by themselves imply ill-posedness. The ratio requires compatibility at characteristic boundaries: in this particular body objective its numerator and denominator contain a cancelling normal factor. One must not infer blow-up from the denominator alone.
What the other cases demonstrate
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.
Quasi-1-D Euler. The speeds explain the two inlet/one outlet primal conditions and reversed dual counts for the stated positive subsonic branch. Counts are necessary, not sufficient: propagation and the boundary constraints together must yield an invertible boundary-value operator. At a sonic point the flux Jacobian is singular, so the current representation using fails; throat compatibility/regularity analysis is required. This does not imply that all transonic adjoints fail to exist. Shocks require shock-position variations and internal conditions. These phenomena are outside the current numerical evidence; see the continuous nozzle analysis of Giles and Pierce. The distinct effects of strong/weak boundary treatment on numerical duals are discussed by Duivesteijn et al..
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.
3. One-dimensional Burgers: second-order reconstruction, continuous adjoints and boundary consistency
Smooth positive branch and the continuous tangent
Consider the steady equation on the fixed interval [0,1], with positive inlet data. Integrating the source from the inlet and choosing the positive square root gives:
As long as the radicand stays strictly positive and bounded away from zero, the smooth solution is unique and depends stably on source and inlet perturbations. Writing v for the state perturbation gives:
The teaching baseline uses s=0.5 and inlet value 1, hence u(x)=sqrt(1+x). The target corresponds to s=sin(πx) and remains fixed during differentiation. Tangent mesh refinement uses the same physical perturbation δs=1, =0.2, rather than a single source interval that shrinks with the mesh.
From variations to the terminal adjoint equation
The objective integrates the squared error against a frozen target. Define the continuous adjoint λ by the following terminal-value problem:
The operator is -uλ′, not -(uλ)′. The tangent differential operator acts on the product uv, and integration by parts directly determines its formal transpose:
Substitution of the tangent equation eliminates v and gives both the continuous total derivative and an explicitly computable adjoint integral:
Thus the continuous sensitivity to a constant source on an interval is the integral of λ over that interval; the inlet sensitivity is λ(0). The pages evaluate these references with composite five-point Gaussian quadrature, independently of the discrete adjoint linear solve.
Well-posedness and discrete inconsistency are different questions
If u≥c>0, the terminal condition uniquely determines λ, with the stability estimate:
The continuous adjoint of the present baseline is therefore well posed. An inlet boundary layer in its discrete counterpart must not be interpreted as ill-posedness of this terminal-value problem. If u approaches or crosses zero, the bound is no longer uniform. Shocks require weak solutions, shock-location variations and corresponding adjoint transmission conditions; this smooth derivation cannot simply be reused.
Second-order reconstruction does not automatically give a globally second-order adjoint
Cell unknowns are cell averages; node unknowns are point values. Second order uses centred interior slopes and one-sided endpoint slopes, reconstructing to the actual face locations before evaluating fluxes. A cell outlet must reconstruct to x=1. The first internal node-based dual face is at h/2, and even its upwind flux must use a trace reconstructed to that location.
Let and denote reconstruction and B denote oriented face scatter. With boundary constants and constraint replacement treated separately, the exact second-order derivatives are:
The discrete Lagrangian is J+ψTR, giving A₂Tψ=-. Balance rows therefore compare the continuous field λ with -ψ. Node row 0 is a strong constraint; ψ₀ is not a control-volume field value. Instead, compare the inlet sensitivities -ψ₀ and λ(0).
The inlet constraint removes its row’s source contribution, while transposed reconstruction propagates weights from other residual rows to near-boundary unknowns. Interior fields, global fields, inlet derivatives and first-source derivatives must therefore be checked separately. Current calculations show nearly second-order interior adjoint errors for second-order reconstruction with the default dissipative inlet, while finite-amplitude local boundary layers reduce global convergence toward order one-half. An upwind cell inlet improves global convergence; the second-order node-based upwind comparison retains a boundary layer.
The first source interval shrinks with the mesh, so its derivative itself tends to zero. A vanishing absolute difference is insufficient evidence of local consistency. The pages additionally report the discrete derivative divided by the continuous interval integral. Total directional-derivative errors are reported separately for a fixed smooth physical perturbation.
GS and GMRES change the solver only
Let P=A₁(u*) be the first-order approximate Jacobian evaluated at the converged second-order state. A fixed number of zero-initialised GS sweeps defines a linear approximate inverse . Tangent and adjoint fixed-point iterations use the exact residual; GMRES also retains the exact second-order operator:
One forward GS sweep gives B₁=(diag(P)+L)⁻¹; the adjoint uses the corresponding reverse sweep B₁T=(diag(P)+LT)⁻¹ on the PT system. Multiple sweeps also use a fixed count and zero initial values so that standard GMRES sees fixed linear preconditioning. Switching to GMRES, or converging linear residuals more tightly, cannot remove adjoint inconsistency caused by discrete boundary closure.
Interactive examples, complete matrices and measured mesh-convergence data are available on the cell-centred page and the node-centred page.
A.2 从方程到可运行参考实现
A.2 From equations to a runnable reference
给定正值 u 和目标 utgt → 分部积分确定 λ(1)=0 → 从 x 到 1 求积 → 比较 −ψ 与源项/入口梯度
Specify positive u and target utgt → derive λ(1)=0 → integrate from x to 1 → compare −ψ and source/inlet gradients
下方代码可以完整复制到本地文件运行。它独立计算连续参考,不调用本页有限体积 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.
node reference_burgers.js展开完整参考代码Show complete reference code
'use strict';
const exact=(which,x,uin=1)=>Math.sqrt(uin*uin+(which==='target'?2*(1-Math.cos(Math.PI*x))/Math.PI:x));
// Composite five-point Gaussian quadrature, also used for exact cell averages.
const gn=[0,.5384693101056831,-.5384693101056831,.906179845938664,-.906179845938664], gw=[.5688888888888889,.4786286704993665,.4786286704993665,.2369268850561891,.2369268850561891];
function integrate(f,a,b,panels=8){let s=0;for(let k=0;k<panels;k++){const h=(b-a)/panels,c=a+(k+.5)*h;for(let j=0;j<5;j++)s+=h/2*gw[j]*f(c+h/2*gn[j]);}return s;}
const lambda=x=>integrate(t=>(exact('base',t)-exact('target',t))/exact('base',t),x,1);
console.log('x,lambda');
for(let i=0;i<=100;i++)console.log(i/100+','+lambda(i/100));
连续伴随、一致性与适定性
Continuous adjoint, consistency and well-posedness
从分部积分推导连续伴随
Derive the continuous adjoint by integration by parts
离散伴随先离散再转置;连续伴随先对微分方程求变分,再分部积分。两条路径不自动等价,尤其在边界附近。这里可以直接手算:对 求变分,得到 。乘以待定函数 后积分,把导数从 转到 。
A discrete adjoint discretises first and transposes afterwards. A continuous adjoint varies the differential equation and then integrates by parts. These routes need not agree automatically, especially at boundaries. Here the calculation is explicit: varying gives . Multiply by an unknown function and integrate to transfer the derivative from to .
为了让最后的体积分等于 ,取 。原始出口的 没有限制,所以其边界系数必须为零,即 。入口扰动是已知设计变化,剩下的边界项给出入口灵敏度。这解释了为什么伴随的已知值在原始出口。
Choose so the last volume integral equals . The primal outlet perturbation is unconstrained, so its boundary coefficient must vanish: . The inlet perturbation is known design data, and its remaining boundary term yields inlet sensitivity. This explains why adjoint data belong at the primal outlet.
“一致性”和“适定性”分别检查什么?适定性问连续问题是否有解、解是否唯一、是否连续依赖数据。一致性比较离散方程及其解能否逼近指定连续问题。局部边界层可以来自离散边界处理;即使每个有限矩阵可逆,也不能仅据此断言连续极限稳定。
Well-posedness asks whether the continuous problem has a unique solution depending continuously on its data. Consistency compares a discretisation with the specified continuous problem. A local boundary layer can come from discrete boundary treatment; invertibility of every finite matrix alone does not establish stability of the continuous limit.
连续 tangent 为 ,。从 出发,选择连续伴随满足下式。积分分部只产生 边界项,不能误写为 。
The continuous tangent satisfies , . Starting from , choose the following continuous adjoint. Integration by parts produces the boundary term ; the adjoint must not be incorrectly written as .
图解:原始信息向右,伴随信息向左
The picture: primal information moves right, adjoint information moves left
原始扰动沿特征速度 从入口向右传播,所以只能在入口给定 ;出口值 是求解结果。伴随方程 在同一特征上反向积分:它从出口出发,所以唯一需要的数据是出口值。边界项 把两者联系起来:出口项含任意的 ,其系数必须为零;入口项含已知的设计扰动,于是成为入口灵敏度。
The primal perturbation travels with the characteristic speed from the inlet to the right, so it can only be prescribed at the inlet, ; the outlet value is part of the solution. The adjoint equation integrates along the same characteristics in reverse, starting at the outlet, so its only datum is an outlet value. The boundary term links the two: the outlet term contains the arbitrary and its coefficient must vanish; the inlet term contains the known design change and becomes the inlet sensitivity.
本例的解、目标与连续伴随
The states and the continuous adjoint of this example
基准源项 、入口值 1 给出 ;目标由 生成,。积分 在出口为零,在入口为 。右图把本页 N=16 的离散伴随画成 (页面的拉格朗日函数是 ,平衡行上 ):内部两种阶数都贴近 ,入口附近的几个点离开曲线,这就是下文分析的边界层。
The base source with inlet value 1 gives ; the target comes from , . The integral vanishes at the outlet and equals at the inlet. The right panel draws this page's N=16 discrete adjoint as (the page's Lagrangian is , so on balance rows): both orders follow in the interior, and a few points near the inlet leave the curve. That departure is the boundary layer analysed below.
逐步推导:tangent 解、伴随积分与对偶恒等式
Step by step: the tangent solution, the adjoint integral and the duality identity
第一步,解 tangent。把 从 0 积分到 ,入口值 直接进入:
Step 1, the tangent. Integrate from 0 to ; the inlet value enters directly:
第二步,解伴随。记 ,则 ;从出口 反向积分:
Step 2, the adjoint. With the equation reads ; integrate backwards from the outlet value :
第三步,不用分部积分也能验证对偶。把 tangent 代入 ,再交换二重积分的次序(积分区域是三角形 ):
Step 3 verifies the duality without integration by parts. Insert the tangent into and exchange the order of the double integral over the triangle :
这正是 :分部积分与交换积分次序是同一个恒等式的两种写法。
This is exactly : integration by parts and the exchange of integration order are two forms of one identity.
第四步,代入本页的方向 、。交换次序后 ,所以伴随给出 ;直接用 积分得到同一个数 。
Step 4 uses this page's direction , . After the exchange, , so the adjoint gives ; integrating directly with gives the same number, .
离散问题有同样的恒等式,而且是代数恒等式:离散设计种子 取同一方向(、),、 给出 。下表对本页格式求出这个方向导数 ,并与连续值比较;同时列出离散 tangent 的最大误差和两种算法之差。
The discrete problem has the same identity, now an algebraic one: with the discrete design seed in the same direction (, ), and give . The table evaluates the directional derivative with this page's scheme and compares it with the continuous value; it also lists the maximum error of the discrete tangent and the difference between the two algorithms.
| 阶数Order | N | ψTRDḊ | 与连续值之差Error | 实测阶Rate | tangent 最大误差Tangent max error | 实测阶Rate | |gTv − ψTRDḊ| |
|---|---|---|---|---|---|---|---|
| 一阶first order | 32 | −0.04270107 | −5.31e−3 | — | 1.04e−2 | — | 2.78e−17 |
| 128 | −0.03889265 | −1.50e−3 | 0.912 | 2.64e−3 | 0.989 | 3.47e−17 | |
| 512 | −0.03777899 | −3.86e−4 | 0.979 | 6.62e−4 | 0.997 | 4.16e−17 | |
| 2048 | −0.03749018 | −9.72e−5 | 0.995 | 1.66e−4 | 0.999 | 0.00e+0 | |
| 二阶second order | 32 | −0.03766552 | −2.73e−4 | — | 3.45e−4 | — | 2.08e−17 |
| 128 | −0.03740981 | −1.68e−5 | 2.009 | 2.35e−5 | 1.938 | 4.16e−17 | |
| 512 | −0.03739403 | −1.05e−6 | 2.003 | 1.51e−6 | 1.980 | 1.04e−16 | |
| 2048 | −0.03739305 | −6.53e−8 | 2.001 | 9.52e−8 | 1.994 | 3.47e−16 |
两种阶数的方向导数都按格式的阶收敛,而离散 tangent 在全域(包括入口)同样按格式的阶收敛。下一小节会看到,离散伴随在入口有一个不随加密消失的边界层;由上面的代数恒等式,同一个方向导数也等于 ,而 没有边界层,所以这个边界层不会降低光滑方向导数的收敛阶(它在二阶格式中的净贡献怎样抵消,见下一小节)。
For both orders the directional derivative converges at the order of the scheme, and so does the discrete tangent over the whole domain, inlet included. The next subsections show that the discrete adjoint carries an inlet boundary layer that does not vanish under refinement; by the algebraic identity above the same directional derivative also equals , and has no layer, so that layer cannot lower the convergence order of derivatives in smooth directions (how its net contribution cancels at second order is shown in the next subsection).
连续伴随是否适定
Is the continuous adjoint well posed?
当前 ,终端值问题唯一且稳定,。如果 穿过零、出现激波或特征方向改变,需要重新分析正则性和边界条件。当前数值边界层本身不能证明连续伴随不适定。
Here , so the terminal-value problem is unique and stable, with . If crosses zero, shocks form or characteristic directions change, regularity and boundary conditions must be reconsidered. The numerical boundary layer here does not itself demonstrate an ill-posed continuous adjoint.
格点第 0 行是约束,其乘子不是控制体内的伴随场值,全域场误差因此排除该行;入口导数单独比较 与 。内部误差固定取 ,全域误差则保留其他边界附近的控制体。
Node row 0 is a constraint, whose multiplier is not an adjoint field value in a control volume, so it is excluded from field-error norms. Its inlet derivative is compared separately with . The interior norm uses the fixed region ; the global norm retains the other near-boundary volumes.
边界附近为什么可能不一致
Why discrepancies can remain near the boundary
默认耗散入口下,二阶格式的内部伴随可接近二阶,但全域 误差会因局部有限幅值边界层接近半阶。格心迎风入口明显改善一致性;格点二阶的重构及强行替换边界仍可保留边界层,不能沿用一阶迎风的结论。首源项导数与连续区间积分的比值还会检出局部不一致,即使光滑方向的总导数已经收敛。
With the default dissipative inlet, second-order interior adjoints can approach order two while the global error approaches order one-half because of a finite-amplitude local boundary layer. An upwind cell inlet improves consistency. Reconstruction and strong boundary-row replacement can retain a boundary layer in the second-order node scheme, so conclusions from first-order upwinding do not carry over automatically. The first source derivative divided by its continuous interval integral reveals local inconsistency even when derivatives in smooth directions converge.
入口伴随边界层:全域误差为什么只有半阶
The inlet adjoint boundary layer: why the global error is of order one half
先看一阶格式的内部行。 的第 行来自 的第 列:左面 贡献 ,右面 贡献 ,对角元是二者之差 :
Start with an interior row of the first-order scheme. Row of is column of : the left face contributes , the right face contributes , and the diagonal is their difference :
代入 ,零阶项相消,剩下 : 正好满足连续伴随 。但这是三项递推,除了光滑解还有第二个解。冻结 、令 ,得到
Substituting cancels the zeroth-order terms and leaves , so satisfies the continuous adjoint . But a three-term recurrence has a second solution besides the smooth one. Freezing and setting gives
本页 、入口 ,所以 。由于 ,第二个解离开入口后每格缩小为原来的 0.231,只占几格宽。格点页的第 0 行是约束 ,它在第 1 行转置中的系数为零,所以第 1 行正好是一般行取 的形式:
On this page and the inlet state is , so . Since , the second solution shrinks by the factor 0.231 per cell away from the inlet and occupies only a few cells. On the node page row 0 is the constraint , whose coefficient in the transposed row 1 is zero, so row 1 is exactly the general row with :
连续伴随在入口不为零(本例 ),光滑解满足不了这个条件,只能加上第二个解。于是在入口附近
The continuous adjoint does not vanish at the inlet (here ), so the smooth solution cannot meet this condition and the second solution must be added. Near the inlet
这就是全域误差的半阶:边界层的幅值 不随加密减小,只有宽度随 收窄;常数 。离开入口几格后 已可忽略,所以内部误差仍按格式的阶收敛。
This is the half order of the global error: the layer amplitude does not shrink under refinement; only its width shrinks with . The constant is . A few cells away from the inlet is negligible, so the interior error still converges at the order of the scheme.
两个直接后果。第一个源区间只进入第 1 行(第 0 行没有源项),权重 :,与连续区间积分 之比趋于 。入口导数却仍然收敛:约束行给出 ,而转置的第 0 行是 ,
Two consequences follow. The first source interval enters row 1 only (row 0 carries no source), with weight : , so its ratio to the continuous interval integral tends to . The inlet derivative still converges: the constraint row gives , and the transposed row 0 is , so
格点迎风对照把第一内部面换成物理通量,一阶时第 1 行变成一般行取 ,第二个解的系数只剩 ;二阶时下面的五项递推仍留下边界层。
The node upwind comparison uses the physical flux on the first interior face; at first order row 1 becomes the general row with , leaving the second solution a coefficient of only . At second order the five-term recurrence below still leaves a layer.
二阶时一个面依赖四个未知量,转置的内部行变成五项递推,冻结系数后的特征方程为
At second order a face depends on four unknowns, so the transposed interior row is a five-term recurrence; with frozen coefficients its characteristic equation is
在入口状态下根为 、 和一对复根 (模 )。复根对离开入口振荡衰减,入口面和端点单边斜率同样给它留下 的系数,所以全域误差仍是半阶,只是常数不同(格点约 0.0171);内部误差和入口导数按二阶收敛。根 属于出口一侧,本例出口端的误差只有 。格点二阶的首源项比值趋于约 0.744,由五项递推与第一对偶面的模板共同决定。
At the inlet state the roots are , and a complex pair of modulus . The complex pair decays with oscillation away from the inlet, and the inlet face and the one-sided end slopes again leave it an coefficient, so the global error is still of order one half with a different constant (0.0171 on this page), while the interior error and the inlet derivative converge at second order. The root belongs to the outlet side; here the outlet error is only . The second-order node first-source ratio tends to about 0.744, set by the five-term recurrence together with the first dual-face stencil.
第三个后果:光滑方向的导数为什么不受影响。设源项密度的设计方向是光滑函数 (种子 )。把 代入离散导数,误差分成边界层与求和公式两部分:
A third consequence: why derivatives in smooth directions are unaffected. Let the design direction of the source density be a smooth function (seed ). Substituting into the discrete derivative splits its error into a layer part and a quadrature part:
格点第 0 行是入口约束,不带源项,求和从 开始,所以求和公式比梯形公式少了入口的半个控制体,误差为 。一阶时边界层的和是 (实测 0.0104),贡献 ,与一阶格式本身的误差同阶。二阶时边界层的和趋于 (实测 -0.0174),它给出 ,正好补回约束行缺少的半个控制体。下表按全部场行计算 :它除以 趋于 ,扣除 后按 缩小。
Node row 0 is the inlet constraint and carries no source, so the sum starts at and misses the inlet half volume of the trapezoidal rule; its error is . At first order the layer sums to (measured 0.0104), an contribution of the same order as the first-order scheme error. At second order the layer sum tends to (measured -0.0174); it contributes and restores exactly the half volume the constraint row lacks. The table evaluates over all field rows: divided by it tends to , and after subtracting it shrinks like .
| 二阶Second order N | Σ Wjej | Σ Wjej / h | (Σ Wjej − hλ(0)/2) / h² |
|---|---|---|---|
| 256 | −7.019e−5 | −0.01797 | −0.147 |
| 1024 | −1.713e−5 | −0.01754 | −0.146 |
| 4096 | −4.255e−6 | −0.01743 | −0.146 |
这种抵消不是巧合。由 ,同一个方向导数也等于 ,而离散 tangent 没有边界层,按格式的阶收敛,所以边界层的净贡献只能按同样的阶消失。只落在第一个区间的方向(首源项)不是光滑方向,这种抵消不会发生,比值仍趋于上面给出的极限。
The cancellation is no accident. Because , the same directional derivative equals , and the discrete tangent has no layer and converges at the order of the scheme; the net contribution of the adjoint layer must therefore vanish at the same order. A direction confined to the first interval (the first source) is not a smooth direction, no such cancellation occurs, and its ratio keeps the limit given above.
下表由本页的离散格式在构建时直接计算(Newton 加带状直接求解,伴随解转置系统)。一阶耗散入口的全域误差与 之比趋于 0.00825,实测阶趋于 0.5;首源项比值趋于 0.385;入口导数误差按一阶减小。
The table is computed at build time from this page's discretisation (Newton with a banded direct solve; the adjoint solves the transposed system). For the first-order dissipative inlet the global error divided by tends to 0.00825 and the measured order to 0.5; the first-source ratio tends to 0.385; the inlet-derivative error decreases at first order.
| 情形Case | N | 全域误差Global error | 全域误差Global error /√h | 实测阶Order | 内部误差Interior error | 首源项比值First-source ratio | 入口导数误差Inlet-derivative error |
|---|---|---|---|---|---|---|---|
| 一阶,耗散入口first order, dissipative | 256 | 8.11e−4 | 0.01297 | — | 5.37e−4 | 0.3986 | −1.27e−3 |
| 1024 | 3.03e−4 | 0.00969 | 0.710 | 1.36e−4 | 0.3882 | −3.20e−4 | |
| 4096 | 1.35e−4 | 0.00864 | 0.583 | 3.40e−5 | 0.3855 | −8.03e−5 | |
| 16384 | 6.52e−5 | 0.00835 | 0.524 | 8.51e−6 | 0.3848 | −2.01e−5 | |
| 二阶,耗散入口second order, dissipative | 256 | 1.07e−3 | 0.01705 | — | 1.93e−6 | 0.7436 | −1.02e−5 |
| 1024 | 5.33e−4 | 0.01706 | 0.499 | 1.20e−7 | 0.7438 | −6.37e−7 | |
| 4096 | 2.67e−4 | 0.01707 | 0.500 | 7.51e−9 | 0.7438 | −3.98e−8 | |
| 16384 | 1.33e−4 | 0.01707 | 0.500 | 4.70e−10 | 0.7439 | −2.49e−9 |
完整连续推导见 中文版连续伴随文档。下方结果均采用当前边界闭合;可在页首切换对照。
See the English continuous-adjoint derivation. The results below use the currently selected boundary closure; switch the comparison at the top of the page.
附录 B:间断、激波位移与离散伴随
Appendix B: discontinuities, shock motion and discrete adjoints
正文先在光滑解上学习离散伴随。本附录进一步问:解出现间断后,怎样构造守恒求解程序,又怎样对这个程序求导?我们使用独立的无源、非定常 Riemann 问题。两种页面在这里都采用同一个单元平均有限体积实现;格点页的正文离散并没有被悄悄替换。附录控制独立于页顶的稳态 GS/GMRES 控制。
The main chapters introduce discrete adjoints using smooth solutions. Here we ask how to build and differentiate a conservative solver when discontinuities occur. This independent experiment solves a source-free, unsteady Riemann problem. Both pages use the same cell-average finite-volume implementation in this appendix; the node-centred discretisation in the main text is unchanged. These controls are independent of the steady GS/GMRES controls above.
B.1 从强形式到熵解:为什么不能直接跨过间断求导
B.1 From the strong form to entropy solutions
在跳跃处,普通意义的空间导数不存在。对一个小控制体积分,再让它包围的激波宽度趋于零,左右通量差必须等于跳跃随界面移动带来的质量变化。因此得到 Rankine–Hugoniot 条件。它给出移动速度,却不能单独选出物理解。Burgers 的熵条件要求压缩激波两侧特征线进入激波;若左值小于右值,物理解是连续展开的稀疏波,而不是满足相同守恒条件的“膨胀激波”。
Ordinary spatial derivatives do not exist at a jump. Integrate over a small control volume enclosing the moving shock: the flux difference must balance the mass change caused by its motion. This gives the Rankine–Hugoniot condition. Conservation determines a speed but does not select the physical solution. Characteristics enter an admissible compressive Burgers shock; increasing Riemann data instead generate a continuous rarefaction, not an expansion shock.
可对照 Clawpack 的 Burgers Riemann 问题。下方包括移动激波 (1,0)、驻定激波 (1,−1) 和稀疏波 (−0.5,0.5)。时间和位置范围保证解析波不碰到边界;左右边界分别以给定的常值外侧状态计算通量。
See the Clawpack Burgers Riemann chapter. The experiment includes a moving shock (1,0), a stationary shock (1,−1), and a rarefaction (−0.5,0.5). The permitted time and position ranges keep the exact waves away from boundaries. Boundary fluxes use the prescribed constant exterior states.
逐步推导:Rankine–Hugoniot 条件
Step by step: the Rankine–Hugoniot condition
守恒律的积分形式对任何控制体都成立,不需要 可微:
The integral form of the conservation law holds for any control volume and does not need to be differentiable:
取 ,两侧分别是常值 、。积分等于两段长度乘各自的值,只有激波位置随时间变化:
Take with constant states and on the two sides. The integral is the two lengths times their values, and only the shock position depends on time:
令两式相等,得到激波速度;对 ,分子可以分解,速度就是两侧状态的平均:
Equating the two expressions gives the shock speed; for the numerator factorises and the speed is the mean of the two states:
熵条件看特征速度 :两侧特征线都要进入激波,即 。代入 ,它等价于 。若 ,同样满足守恒的“膨胀激波”两侧特征线离开激波,不是物理解,物理解是下图 (c) 的稀疏波扇区。
The entropy condition concerns the characteristic speed : characteristics from both sides must enter the shock, . With this is equivalent to . For an "expansion shock" would conserve mass but its characteristics leave it; the physical solution is the rarefaction fan of panel (c) below.
B.2 网格、初值投影与限幅重构
B.2 Mesh, initial projection and limited reconstruction
未知量是每个单元的平均值 ui,初始时刻记作 ui⁰。若初始跳跃切过一个单元,就按左右部分长度加权,不能简单按单元中心的位置赋值,否则小幅改变 a 时整个离散初值不变,计算出的导数会错误地归零。下面 θ 是左状态在该单元所占的比例。
The unknowns are the cell averages ui, written ui⁰ at the initial time. A cell cut by the initial jump receives a length-weighted average. Sampling the state only at the cell centre would make sufficiently small changes in a invisible, giving a misleading zero derivative. Here θ is the fraction occupied by the left state.
当 a 正好在网格面上,上式不再给出唯一的双侧导数:a 向左或向右移动时变化的是不同单元。实现选取右侧方向作为初始 tangent,界面会明确提示,并同时列出单侧差分。离散初值的向量导数随网格加密变得更尖,并不应期待它逐点收敛为普通函数。
At a mesh face there is no unique two-sided derivative of this projection: moving a to the left or to the right changes different cells. The implementation selects the right-side initial direction, displays a warning, and reports both one-sided differences. The initial tangent becomes more concentrated on refinement; pointwise convergence to an ordinary function is not expected.
一阶令 d=0。限幅模式在两侧差分异号或有零时令 d=0,否则取绝对值较小的差分;相等时选左差分。首末单元不重构。它在光滑、非极值区域允许二阶重构,在间断附近自动降低精度;不能把“二阶模式”理解成含激波解在所有范数下都有二阶收敛。下面的图保留每个网格面和每个单元平均值,选中单元还能查看具体限幅运算。
First order sets d=0. The limited mode chooses zero for opposite signs or a zero difference; otherwise it selects the smaller magnitude, choosing the left difference at a tie. Endpoint cells are not reconstructed. Second-order reconstruction is possible in smooth non-extremal regions, while accuracy is reduced near a discontinuity. This does not imply second-order convergence in every norm for solutions with shocks. The plot shows every mesh face and cell average; selecting a cell exposes its limiter calculation.
逐步推导:切割单元的初值导数
Step by step: the derivative of the cut-cell average
设跳跃 落在第 个单元 内。单元平均把两段分开积分:
Let the jump lie in cell , . The cell average integrates the two pieces separately:
其他单元不含 ,导数为零。于是 :离散导数是一个高 、宽 的尖峰,总质量与连续导数 相同。
No other cell contains , so their derivatives vanish. Hence : the discrete derivative is a spike of height and width whose total mass equals that of the continuous derivative .
B.3 守恒通量与实际时间推进
B.3 Conservative flux and implemented time integration
这是局部 Lax–Friedrichs 形式的数值通量(和正文一样记作 H,两侧迹记作 u⁻、u⁺),本实验把波速上界 α 固定为 1.2,而非每一步计算 max|u|。给定算例的状态范围不超过 1,因此这个选择覆盖波速,也避免把 max/绝对值的分支变化混入求导。每个内部面只计算一次通量,并以相反符号进入相邻单元;所以质量变化精确等于边界净通量,只有舍入误差。
This Lax–Friedrichs-form flux (written H as in the main text, with traces u⁻ and u⁺) uses a fixed wave-speed bound α=1.2 rather than recomputing max|u|. The supplied states have magnitude at most one, so this covers their wave speeds and avoids differentiating another max/absolute-value branch. Each internal face contributes with opposite signs to adjacent cells, making mass change equal the net boundary flux up to roundoff.
使用两级 SSP-RK2;步数取整后重新设置 Δt,恰好到达终止时间 T。改变设计参数 a 不改变步数、Δt、网格或 α。这样比较的是同一个离散程序。这里是显式非定常算法,没有需要 GS 求解的稳态线性方程;若改为隐式推进,每个时间层仍可用正文的 GS/预处理 GMRES 求解,但还必须把所有时间层的依赖反向传播,不能只求一个最终稳态伴随。
The implementation uses two-stage SSP-RK2. An integer number of steps is chosen and Δt adjusted to reach T exactly. Changing a leaves the mesh, time schedule and α unchanged, so comparisons differentiate the same discrete program. This explicit transient algorithm has no steady linear system to solve with GS. An implicit version could use GS/preconditioned GMRES at each time level, but its adjoint must still reverse dependencies through all time levels.
这里 N 表示单元数,K 表示时间步数;本附录的 di 是单元内增量(等于斜率乘 h),与正文按长度求导的斜率符号不同。以下示意图标出两条必须同时保留的 RK 依赖路径。
Here N counts cells and K counts time steps. The appendix uses di for an increment (slope times h), unlike the spatial derivative used for slopes in the main text. The diagram shows both RK dependency paths that must be retained.
B.4 Matrix-free tangent:求导到每个选中的分支
B.4 Matrix-free tangent through each selected branch
先保存每个 RK 子阶段选中的 minmod 分支。左分支给出 ,右分支给出 ,零分支给出 。再求左右重构值、面通量和单元残差的切向量。这是局部运算得到的 ,不组装全局矩阵。注意:保存分支仍要对该分支里的差分求导;把整个限幅斜率当常数冻结通常不是原离散程序的导数。
Save the selected minmod branch at each RK stage. Its slope tangent is for the left branch, for the right branch, and for the zero branch. Carry the tangents through the reconstructed traces, face fluxes and residual scatter. These local operations evaluate without a global matrix. Saving a branch still requires differentiating the difference inside it; freezing the whole limited slope generally differentiates a different program.
在 minmod 的零点或等幅交界处,所保存的是一个选定分支的线性映射,不保证存在对所有方向都适用的经典导数。即使 tangent 与 adjoint 的点积完全相等,也只证明两段代码互为转置;还需要比较单侧差分、分支切换以及步长变化。不能直接把正文的复数步长算法套到比较运算和 minmod 上。
At zero or equal-magnitude limiter boundaries, the saved map is a selected-branch linearisation, not necessarily a classical derivative valid for every direction. A perfect tangent/adjoint dot test only establishes transposition of those maps. One-sided differences, branch changes and step-size dependence still matter. The complex-step procedure in the main text cannot be applied blindly through comparisons and minmod.
B.5 反向时间 adjoint:转置面循环和 RK 计算图
B.5 Reverse-time adjoint: face loops and the RK graph
对一个面,从空间算子的种子 出发,面通量的伴随是 ;边界外侧种子为零。乘以通量对 、 的偏导,累加到两侧单元的伴随,并把 的贡献传给斜率的伴随 。最后按照保存的 minmod 分支,把 加到对应的两个单元。这直接计算 。随后的 RK 反传必须同时保留 到最终平均值的直接路径与经过 的路径。
For each face, start from the seed of the spatial operator: the face-flux adjoint is , with zero exterior seeds. Multiply by the flux derivatives with respect to and , accumulate into the adjoints of the two adjacent cells, and pass the contributions to the slope adjoints . Reverse the saved slope branch into its two cells. This evaluates directly. Reversing RK must include both the direct averaging path from and the path through .
本实验存储每个时间层和两级重构分支,内存随“单元数×时间步数”增长。工业非定常伴随通常需要检查点和重算。图中画的是 ū/h,便于与连续伴随场比较;真正用于点积的是含 h 的 ū,不能再额外乘一次 h。
This experiment stores states and both reconstruction stages, using memory proportional to cells times steps. Large unsteady adjoints usually need checkpointing and recomputation. The plot displays ū/h for comparison with a continuous adjoint density; dot products use ū itself, which already contains h.
u = exact_cell_average(step(a))
for n = 0,...,K-1:
A = stage(u) # flux + saved limiter branch
y = u + dt*A.L
B = stage(y)
u = (u + y + dt*B.L)/2
save(A, B)
u_b = h*q*Gprime(u) # u_b = ū^K: the seed dJ_h/du^K
for n = K-1,...,0: # reverse order: the adjoint of y first
y_b = (u_b + dt*vjp(B[n],u_b))/2
u_b = u_b/2 + y_b + dt*vjp(A[n],y_b)
gradient = dot(u_b, initial_direction(a))逐步推导:一个面的反向传播
Step by step: reversing one face
前向计算中,内部面 的通量 只进入两个单元的 。以 为种子,它的伴随收集这两条路径:
In the forward computation, the flux of the interior face enters of only two cells. With the seed , its adjoint collects these two paths:
再乘以通量对两侧迹的偏导,然后按重构 、 的转置累加(反向语句按与前向相反的顺序执行):
Multiply by the flux derivatives with respect to the two traces, then accumulate through the transposed reconstruction , (the reverse statements run in the opposite order of the forward ones):
最后把每个 按保存的分支送回两个单元:左分支 给出 、;右分支类推;零分支什么也不做。边界面的外侧状态是常数,没有对应的伴随。
Finally send each back to two cells through the saved branch: the left branch gives and ; the right branch is analogous; the zero branch does nothing. Exterior states at boundary faces are constants and have no adjoint.
逐步推导:SSP-RK2 的反向与点积守恒
Step by step: reversing SSP-RK2 and the conserved dot product
一步 SSP-RK2 的 tangent 是两个线性映射的组合:
The tangent of one SSP-RK2 step combines two linear maps:
反向时,横线量是目标对该中间量的导数(伴随),语句按相反顺序执行:先算 ,再算 。 只通过第二式影响 ; 既直接进入 ,又通过 进入,所以两条路径相加:
In reverse, a barred quantity is the derivative of the objective with respect to an intermediate quantity (its adjoint), and the statements run in the opposite order: first , then . affects only through the second map; enters both directly and through , so the two paths add:
用这两式可以直接证明每一步的点积不变(AD 的对偶恒等式),于是 ,即 adjoint 与 tangent 给出同一个 :
These two formulas prove directly that each step preserves the dot product (the duality identity of AD), so : the adjoint and the tangent give the same :
B.6 激波位移导数与连续伴随的适定性
B.6 Shock displacement and continuous-adjoint well-posedness
改变初始位置 a 会移动激波。把解写成两个常值和 Heaviside 阶跃,就会得到 Dirac 分布,而不是一个有界的普通 tangent 场。目标函数仍可有有限导数:应先对移动界面的积分求导,不能在激波处随意代入一个 G′(u)。以下方括号表示左值减右值。
Changing a moves the shock. Differentiating the Heaviside representation gives a Dirac distribution rather than a bounded tangent field. An integral objective can nevertheless have a finite derivative: differentiate the moving-interface integral instead of assigning an arbitrary G′(u) at the jump. Brackets below denote left minus right.
这个推导也解释了为何不能写成 G′(u(xs))·[u]:间断点的 u 没有唯一值,非线性 G 的跳跃是 [G]。对二次目标,。若忽略激波位移,即使两侧光滑区的线性化都正确,也会漏掉整个位置导数。
This derivation explains why G′(u(xs))·[u] is not a valid replacement: u has no unique value at the jump, and the nonlinear jump is [G]. For a quadratic objective, . Ignoring shock motion loses the entire position derivative even when both smooth-side linearisations are correct.
在两侧光滑区域形式上有 −λt−uλx=0,终端值为 qG′。但这不足以定义跨激波的伴随:必须包含界面移动的贡献。对此常值 Riemann 激波,把进入激波的初始特征线区域识别出来,得到下面的初始伴随平台。它满足 λ(a,0)·[u]=q(xs)[G],从而恢复位置导数。交互在初始时刻画出这个解析参考;稀疏波没有激波平台,不画此参考。
In smooth regions the formal adjoint is −λt−uλx=0 with terminal data qG′. These equations alone do not define an adjoint across the shock: interface motion must be included. For this constant-state Riemann shock, identifying initial characteristics captured by the shock gives the plateau below. Its value times [u] recovers q(xs)[G]. The interaction overlays this analytic reference at the initial time; no shock plateau is drawn for a rarefaction.
稀疏波可直接在扇区内求导:dJ/da=−∫扇区 q(x)G′((x−a)/T)/T dx,边缘因状态连续而没有额外跳跃项。驻定激波 (1,−1) 的二次目标有 [G]=0,所以连续位置导数恰为零;数值上涂抹的激波可能产生非零的位置灵敏度。这里的差异是需要研究的结果,不能通过修改解析答案来“对齐”。
For a rarefaction, dJ/da=−∫fan q(x)G′((x−a)/T)/T dx; continuity at fan edges eliminates jump contributions. For the stationary (1,−1) shock and quadratic objective, [G]=0 and the continuous position derivative is exactly zero. A numerically smeared shock may produce nonzero position sensitivity. Such discrepancies are a result to investigate, not a reason to adjust the analytic reference.
逐步推导:移动界面的求导(Leibniz 公式)
Step by step: differentiating a moving interface (Leibniz rule)
激波情形的目标由两段常值积分组成,积分上下限随 移动。对积分限求导(Leibniz 公式),被积函数本身与 无关:
For a shock the objective consists of two integrals of constants whose limits move with . Differentiating the limits (Leibniz rule), with integrands independent of :
本页默认 、:移动激波的 , 时 , 时为 ;驻定激波 , 分别为 2 与 0,得到 与 0。
With this page's defaults , the moving shock has , giving for and for ; the stationary shock has and equal to 2 and 0, giving and 0.
逐步推导:伴随平台来自质量守恒
Step by step: the adjoint plateau from mass conservation
伴随值就是目标对初值局部扰动的灵敏度。在初始时刻 处加一小块质量 ,即 。它沿特征线 平移。若在 之前没有碰到激波,终端的目标变化是 ,这就是平台两侧的值。若它先进入激波,扰动被激波吸收,两侧状态恢复为 、,多出的质量只能表现为激波位置的移动:
An adjoint value is the sensitivity of the objective to a local change of the initial data. Add a small mass at at the initial time, . It travels along the characteristic . If it does not meet the shock before , the terminal change of the objective is : these are the values on the two sides of the plateau. If it enters the shock first, the shock absorbs it, both states return to and , and the extra mass can only appear as a shift of the shock:
这个值与 无关,只要 的特征线在 之前进入激波。左侧特征线 在 恰好到达激波的起点是 ,右侧是 ,两者之间就是平台。 时平台与两侧连续(值 1.446); 时平台值 0.723,两侧分别是 与 0,出现两个跳跃。平台值乘 正好是上面的 。
The value does not depend on , provided the characteristic from enters the shock before . The left characteristic reaches the shock exactly at from , the right one from ; the plateau lies between them. For the plateau is continuous with both sides (value 1.446); for the plateau value is 0.723 while the sides are and 0, with two jumps. The plateau value times is exactly the above.
逐步推导:稀疏波的导数
Step by step: the derivative of the rarefaction
稀疏波有三段:扇区外两段常值,扇区内 。三段的积分限都随 移动:
The rarefaction has three pieces: constants outside the fan and inside it. All the limits move with :
扇区边缘状态连续,,所以 Leibniz 公式产生的四个边界项两两抵消,只剩被积函数对 的导数;再令 :
The state is continuous at the fan edges, , so the four boundary terms of the Leibniz rule cancel in pairs and only the derivative of the integrand remains; then set :
对 、: 时为 ; 时奇次项相消,只剩 。这些值就是 B.7 实验结果表中“连续位置导数”一行的数。
For and : gives ; for the odd terms cancel and only remains. These are the values in the “Continuous position derivative” row of the B.7 results table.
还要区分另一种适定性问题:对稳态、无源的 (1,−1) 问题,左右通量相同,很多激波位置都满足相同边界数据,位置没有被唯一选定。其平移方向可能对应零模,稳态线性化可能奇异。这里的非定常问题通过初始位置 a 选定轨迹,因此不能用稳态的不唯一性否定当前离散时间程序的求导。
A separate well-posedness issue arises for the source-free steady (1,−1) problem. Equal boundary fluxes admit multiple shock positions; translation may form a null mode and make the steady linearisation singular. The transient initial position a selects a trajectory. Steady nonuniqueness therefore does not invalidate differentiation of this specified discrete time program.
Giles 与 Ulbrich(2010,第二部分)研究了激波处离散伴随的收敛,并展示了原始解收敛而伴随平台不正确的情形;其证明针对特定数值耗散及随加密增长的激波覆盖点数。本附录的固定 α、minmod 方法不自动满足该定理的假设。因此下方同时显示离散验证和连续参考误差,不把前者当作后者的证明。
Giles and Ulbrich (2010, Part II) analyse discrete-adjoint convergence at shocks and exhibit incorrect adjoint plateaux despite primal convergence. Their proof uses particular dissipation scaling and an increasing number of points across the shock under refinement. That theorem does not automatically cover this fixed-α, minmod scheme. The experiment reports discrete verification separately from error against the continuous reference.
B.7 可执行实验:网格、限幅、时间反传和差分
B.7 Executable experiment: mesh, limiter, reverse time and FD
误差表中的 L1 是 hΣ|ui−精确单元平均i|,不是对分片常值数值解与点态精确解之差的连续积分;实测阶为 log₂(/)。两种误差在切过激波的单元中不同。图中的连续伴随参考按单元中心采样。分支改变数统计的是任一 RK 子阶段发生变化的时间层,而不是改变的单元总数。
The tabulated L1 error is hΣ|ui−exact cell averagei|, not the continuous integral of the error between a piecewise-constant numerical field and the pointwise exact solution. Observed order is log₂(/); these error definitions differ in shock-cut cells. The continuous-adjoint reference is sampled at cell centres. Branch-change counts refer to time levels with a change in either RK stage, not the total number of switched cells.
使用顺序:先比较两种空间离散的原始解和网格细化误差;再拖动时间层,看 tangent 如何集中到激波附近、伴随如何反向形成平台;最后改变 a(也可输入恰好位于网格面的数),比较单侧差分与中心差分。限幅分支改变的时间层数仅作诊断:小数值尾部的分支变化不一定显著影响目标,但不能忽略真实的不光滑点。
Suggested sequence: compare primal solutions and refinement errors for both spatial schemes; move the time slider to see tangent concentration and reverse-time adjoint evolution; then vary a, including positions exactly on mesh faces, and compare one-sided and central differences. The count of time levels with changed limiter branches is diagnostic: switches in tiny numerical tails need not affect the objective significantly, but actual nonsmooth points cannot be ignored.
B.8 小结与实现范围
B.8 Summary and implementation scope
本附录完整实现了三个 Riemann 问题的守恒时间推进、初值位置 tangent、反向时间离散 adjoint,以及解析解、质量平衡、点积和步长扫描验证。它没有实现任意激波相互作用的连续伴随,也没有证明当前耗散下伴随随网格加密收敛。教学的主线仍是对实际离散程序建立可核查的伴随;简单连续解用于验证和揭示局限。走向工业非结构网格 RANS 时,同样要交代限幅、边界、湍流闭合和求解过程的具体求导约定,不能假定存在一个可直接作为严格基准的完整连续伴随。
This appendix implements conservative time integration, initial-position tangents, reverse-time discrete adjoints, and analytic-solution, mass-balance, dot-product and step-size checks for three Riemann problems. It does not implement continuous adjoints for arbitrary shock interactions or prove adjoint convergence for this dissipation. The teaching focus remains an auditable adjoint of the actual discrete program. Simple continuous solutions expose both agreement and limitations. Industrial unstructured RANS additionally requires explicit differentiation conventions for limiters, boundaries, turbulence closure and solver operations; a complete continuous adjoint cannot simply be assumed to provide a rigorous reference.
查看本实验完整数值实现(JavaScript)Full numerical implementation (JavaScript)
/* Appendix B: conservative unsteady cell averages, selected-branch JVP/VJP. */
(function(root){'use strict';
const zeros=n=>Array(n).fill(0),dot=(a,b)=>a.reduce((s,x,i)=>s+x*b[i],0);
const cases={shock:[1,0],stationary:[1,-1],fan:[-.5,.5]};
function model({n=80,T=.15,a=.371,kind='shock',order=1,alpha=1.2,objective='linear'}={}){
if(![n,T,a,alpha].every(Number.isFinite)||!['linear','quadratic'].includes(objective)||!cases[kind]||![1,2].includes(order)||!Number.isInteger(n)||n<8||T<0||T>.25||a<.2||a>.6||alpha<1)throw Error('Invalid shock experiment');
const [left,right]=cases[kind],h=1/n,x=Array.from({length:n},(_,i)=>(i+.5)*h),steps=T?Math.ceil(T/(.35*h/alpha)):0,dt=steps?T/steps:0;
const G=u=>objective==='linear'?u:.5*u*u,Gp=u=>objective==='linear'?1:u,q=x=>1+x;
function coordinate(position){const z=position/h,r=Math.round(z);return Math.abs(z-r)<=8*Number.EPSILON*Math.max(1,Math.abs(z))?r:z;}
function initial(position=a){const z=coordinate(position);return x.map((_,i)=>right+(left-right)*Math.max(0,Math.min(1,z-i)));}
function initialDirection(position=a){const v=zeros(n);v[Math.min(n-1,Math.floor(coordinate(position)))]=(left-right)/h;return v;}
function stage(u){
const slope=zeros(n),weights=Array.from({length:n},()=>[]),branch=Array(n).fill('zero');let ties=0;
if(order===2)for(let i=1;i<n-1;i++){
const l=u[i]-u[i-1],r=u[i+1]-u[i];
if(l*r<=0){if(l===0||r===0)ties++;continue;}
const pick=Math.abs(l)<=Math.abs(r);if(Math.abs(l)===Math.abs(r))ties++;
slope[i]=pick?l:r;weights[i]=pick?[[i,1],[i-1,-1]]:[[i+1,1],[i,-1]];branch[i]=pick?'left':'right';
}
const flux=[],ca=[],cb=[];
for(let f=0;f<=n;f++){const L=f?u[f-1]+.5*slope[f-1]:left,R=f<n?u[f]-.5*slope[f]:right;flux[f]=.25*(L*L+R*R)-alpha/2*(R-L);ca[f]=L/2+alpha/2;cb[f]=R/2-alpha/2;}
return {u,slope,weights,branch,ties,flux,ca,cb,L:u.map((_,i)=>(flux[i]-flux[i+1])/h)};
}
function jvp(t,v){const ds=t.weights.map(w=>w.reduce((s,[j,c])=>s+c*v[j],0)),df=[];for(let f=0;f<=n;f++)df[f]=t.ca[f]*(f?v[f-1]+.5*ds[f-1]:0)+t.cb[f]*(f<n?v[f]-.5*ds[f]:0);return v.map((_,i)=>(df[i]-df[i+1])/h);}
function vjp(t,w){const y=zeros(n),sb=zeros(n);for(let f=0;f<=n;f++){const eta=((f<n?w[f]:0)-(f?w[f-1]:0))/h;if(f){const b=t.ca[f]*eta;y[f-1]+=b;sb[f-1]+=.5*b;}if(f<n){const b=t.cb[f]*eta;y[f]+=b;sb[f]-=.5*b;}}for(let i=0;i<n;i++)for(const [j,c] of t.weights[i])y[j]+=c*sb[i];return y;}
function run(position=a,keep=true){let u=initial(position);const history=keep?[u.slice()]:[],tape=[],branches=[];let conservation=0;
for(let k=0;k<steps;k++){
const A=stage(u),y=u.map((v,i)=>v+dt*A.L[i]),B=stage(y),next=u.map((v,i)=>.5*v+.5*(y[i]+dt*B.L[i]));
const dm=h*next.reduce((s,v,i)=>s+v-u[i],0),boundary=dt/2*(A.flux[0]-A.flux[n]+B.flux[0]-B.flux[n]);conservation=Math.max(conservation,Math.abs(dm-boundary));
branches.push(A.branch.join(',')+'|'+B.branch.join(','));if(keep)tape.push({A,B});u=next;if(keep)history.push(u.slice());
}
const J=h*u.reduce((s,v,i)=>s+q(x[i])*G(v),0);return {u,J,tape,history,branches,conservation};
}
function derivatives(run){let v=initialDirection(),tangent=[v.slice()];
for(const {A,B} of run.tape){const av=jvp(A,v),vy=v.map((z,i)=>z+dt*av[i]),bv=jvp(B,vy);v=v.map((z,i)=>.5*z+.5*(vy[i]+dt*bv[i]));tangent.push(v.slice());}
let w=run.u.map((u,i)=>h*q(x[i])*Gp(u));const dual=Array(steps+1);dual[steps]=w.slice();
for(let k=steps-1;k>=0;k--){const {A,B}=run.tape[k],bw=vjp(B,w),wy=w.map((z,i)=>.5*z+.5*dt*bw[i]),aw=vjp(A,wy);w=w.map((z,i)=>.5*z+wy[i]+dt*aw[i]);dual[k]=w.slice();}
return {tangent,dual,forward:dot(run.u.map((u,i)=>h*q(x[i])*Gp(u)),v),adjoint:dot(w,initialDirection())};
}
const speed=(left+right)/2;
function exact(xx,time=T,position=a){if(time===0)return xx<position?left:right;if(left>right)return xx<position+speed*time?left:right;const z=(xx-position)/time;return z<left?left:z>right?right:z;}
function integral(f,lo,hi){if(hi<=lo)return 0;const c=(lo+hi)/2,r=(hi-lo)/2,z=Math.sqrt(3/5);return r*(5/9*f(c-r*z)+8/9*f(c)+5/9*f(c+r*z));}
function exactIntegral(f,lo,hi,time=T,position=a){const breaks=[lo,hi,...(left>right||time===0?[position+speed*time]:[position+left*time,position+right*time])].filter(v=>v>=lo&&v<=hi).sort((a,b)=>a-b);let s=0;for(let i=1;i<breaks.length;i++)s+=integral(xx=>f(xx,exact(xx,time,position)),breaks[i-1],breaks[i]);return s;}
const exactAverages=time=>x.map((_,i)=>exactIntegral((xx,u)=>u,i*h,(i+1)*h,time)/h);
function referenceGradient(){if(left>right||T===0)return q(a+speed*T)*(G(left)-G(right));return -integral(xx=>q(xx)*Gp((xx-a)/T)/T,a+left*T,a+right*T);}
function continuousInitialDual(xx){if(left<=right)return null;const lo=a+(speed-left)*T,hi=a+(speed-right)*T;if(xx<lo)return q(xx+left*T)*Gp(left);if(xx>hi)return q(xx+right*T)*Gp(right);return q(a+speed*T)*(G(left)-G(right))/(left-right);}
return {n,T,a,kind,order,alpha,objective,h,x,dt,steps,left,right,speed,G,Gp,initial,initialDirection,stage,jvp,vjp,run,derivatives,exact,exactAverages,referenceGradient,continuousInitialDual,exactJ:()=>exactIntegral((xx,u)=>q(xx)*G(u),0,1),aligned:Number.isInteger(coordinate(a))};
}
const api={model,dot};if(typeof module!=='undefined')module.exports=api;else root.BurgersShock=api;
})(typeof window!=='undefined'?window:globalThis);
附录 C:显式矩阵与 LU 直接伴随求解
Appendix C: explicit matrices and direct LU adjoint solves
C.1 一阶和二阶都可以直接求解
C.1 Direct solves for both first and second order
本附录不采用 matrix-free 求解。先按顶端选定的一阶或二阶离散,显式组装完整 Jacobian A,再存储 AT 并进行带部分主元的 LU 分解,最后以前代和回代求出伴随。默认只有 8 个格心未知量或 9 个格点未知量;即使切换到 32 个网格区间,稠密直接法的成本也很小。正文的 matrix-free GS/GMRES 保留用于理解大型问题的求解方式。
This appendix does not use a matrix-free solve. Explicitly assemble the complete Jacobian A for the selected first- or second-order discretisation, store AT, factor it with partial-pivot LU, and compute the adjoint by forward/back substitution. The default grid has only eight cell or nine node unknowns; even 32 intervals remain a small dense problem. The main matrix-free GS/GMRES path remains available to explain methods for larger systems.
C.2 按面显式组装完整矩阵
C.2 Explicit face-by-face matrix assembly
每个面先求左右重构权重和通量偏导 cL、cR,把 cL·左权重+cR·右权重加到左残差行、减到右残差行。二阶权重会延伸到更远的邻点,必须保留全部项。格点入口最后把第 0 行替换为 [1,0,…]。源项不依赖状态,因而不贡献 A;它仍通过 进入最终设计梯度。这里调用完整矩阵组装 matrix(u,false),不使用一阶预处理近似。
For each face, evaluate reconstruction weights and flux partials cL,cR. Add cL times left weights plus cR times right weights to the left residual row and subtract them from the right row. Retain all extended second-order neighbours. On the node grid, replace row zero with [1,0,…] after assembly. The state-independent source contributes no A entries but still enters the design gradient through . Use the complete assembly matrix(u,false), not the first-order preconditioner.
A = zeros(n,n)
for face f:
q = face_values(u,f)
for (row,sign) in [(left,+1),(right,-1)]:
if row exists:
for (j,w) in left_reconstruction: A[row,j] += sign*q.cL*w
for (j,w) in right_reconstruction: A[row,j] += sign*q.cR*w
if node_centred: A[0,:] = [1,0,...]
M = transpose(A)
LU, permutation = factor_with_partial_pivoting(M)
psi = triangular_solve(LU, permutation, -J_u)C.3 LU 分解、行置换、前代与回代
C.3 LU factorisation, permutations and triangular solves
每列从剩余行中选绝对值最大的主元,交换整行,并记录置换 Π。消元乘数保存在严格下三角 L(其对角线为 1),更新后的上三角保存 U。整行交换包括先前已经存储的 L 乘数;右端项也必须采用同一个置换。遇到数值奇异主元时报告失败,不人为替换零主元。不形成逆矩阵。
Choose the largest-magnitude remaining pivot in each column, swap whole rows and record Π. Store elimination multipliers below the diagonal (L has unit diagonal) and U on/above it. Full row swaps include previously stored L multipliers. Apply the same permutation to the right-hand side. Report numerically singular pivots instead of replacing zeros arbitrarily. Do not form a matrix inverse.
下面的网格值、精确矩阵、消元步骤、置换右端、前代结果 y 和伴随 ψ 均现场计算。切换顶端阶数、网格数或入口闭合后自动重新组装和分解。消元步骤不是外层迭代,不存在 GS 意义下的收敛历史。
The mesh values, exact matrix, elimination steps, permuted right-hand side, intermediate y and adjoint ψ are computed live. Changing order, grid size or inlet closure automatically reassembles and refactors. Elimination stages are not outer iterations and do not form a GS convergence history.
逐步推导:一次分解同时服务 tangent 与 adjoint
Step by step: one factorisation serves both the tangent and the adjoint
本页分解的是 。如果改为分解 ,同一组因子也能解伴随:由 得 ,其中 是下三角、 是单位上三角,所以伴随只是把两次三角求解的顺序和方向换过来:
This page factors . Had it factored , the same factors would also solve the adjoint: gives , where is lower triangular and unit upper triangular, so the adjoint merely swaps the order and direction of the two triangular solves:
分解本身约占 次运算,每次三角求解对只需约 次。多个目标函数(多个右端)共用一个 时,只分解一次。
The factorisation costs about operations and each pair of triangular solves only about . Several objectives (several right-hand sides) sharing one need one factorisation.
带状矩阵:本页的消元没有交换行,也没有填充
Banded matrices: on this page elimination swaps no rows and creates no fill
设 的下、上带宽为 、(一阶 ,二阶 )。第 步只有 个乘数 非零,更新只触及 个元素,所以全程约 次运算。部分主元若交换行, 的上带宽最多变为 。对本页的矩阵(格点,一阶与二阶, 与 ),构建时的分解一次也没有交换行, 与 的带宽与 相同,填充为零。这是对这些矩阵的计算结果,不是一般保证。
Let have lower and upper bandwidths and ( at first order, at second order). Step has only nonzero multipliers and updates only entries, so the whole factorisation costs about operations. If partial pivoting swaps rows, the upper bandwidth of can grow to . For this page's matrices (node, first and second order, and ) the build-time factorisation never swapped a row: and keep the bandwidths of and the fill is zero. This is a computed fact about these matrices, not a general guarantee.
C.4 小网格的速度、验证和适用范围
C.4 Small-grid performance, verification and scope
这里直接法求解阶段完全不调用 matrix-free 算子;求解完成后才用正文独立的反向算子检查真残差 ATψ+。同时比较 GS、GMRES 的伴随和全部设计梯度。计时包含直接法的精确矩阵组装、转置、LU 和三角求解;迭代法包含低阶预处理组装及求解。共享已收敛流场,不计原始求解和绘图,每种方法预热后测 7 次取中位数。
The direct solve itself never calls a matrix-free operator. Only afterwards does the main reverse operator independently check ATψ+. Compare adjoints and all gradients with GS and GMRES. Timing includes exact assembly, transpose, LU and triangular solves for the direct method, and low-order setup plus solution for iterative methods. Share the converged flow, exclude primal solves and rendering, warm up each method and report the median of seven runs.
稠密存储为 O(n²)、LU 分解为 O(n³)、复用分解后每个新右端为 O(n²)。因此小系统直接法通常很快,但应以实测结果判断;计时分辨率、浏览器和求解器参数都会影响排名。若改变目标但保持同一个 A,可以复用 LU;若改变离散阶数或状态则必须重分解。工业大规模非结构网格系统需要考虑稀疏填充和内存,不能照搬稠密直接法。小后向误差也不等于病态系统具有小前向误差。
Dense storage costs O(n²), LU O(n³), and each new right-hand side with reused factors O(n²). Direct methods are often fast here, but the live measurements determine the ranking, which depends on timing resolution, browser and solver settings. Reuse LU for a new objective with unchanged A; refactor when order or state changes. Large industrial unstructured systems require attention to sparse fill and memory. Small backward error also does not guarantee small forward error in an ill-conditioned system.
逐步推导:后向误差为什么能换算成前向误差
Step by step: turning the backward error into a forward error
Rigal–Gaches 定理说明,上式的 恰好是使 成为某个扰动系统精确解的最小相对扰动:存在 、,
The Rigal–Gaches theorem says that above is exactly the smallest relative perturbation for which solves a perturbed system: there are and with
与精确解 相减,得 。两边乘 并取范数:
Subtracting the exact equation gives . Multiply by and take norms:
最后一步用了 。本页默认网格(格点,):一阶 、;二阶 、。所以伴随的相对前向误差不超过约 。条件数随网格加密增大( 时分别为 88 与 142),小后向误差并不自动意味着小前向误差。
The last step used . On this page's default mesh (node, ): first order and ; second order and . The relative forward error of the adjoint is therefore at most about . The condition number grows under refinement (88 and 142 at ), so a small backward error does not by itself mean a small forward error.
完整直接法实现Complete direct implementation
/* Explicit Burgers Jacobian assembly and partial-pivot LU: no matrix-free solve. */
(function(root){'use strict';const K=typeof module!=='undefined'?require('./core_1d'):root.FV1;
function factor(M,trace=false){const n=M.length,lu=M.map(r=>r.slice()),permutation=Array.from({length:n},(_,i)=>i),steps=[],scale=Math.max(...M.flat().map(Math.abs));if(!scale||!Number.isFinite(scale))throw Error('Invalid or zero matrix');for(let k=0;k<n;k++){let pivot=k;for(let i=k+1;i<n;i++)if(Math.abs(lu[i][k])>Math.abs(lu[pivot][k]))pivot=i;if(Math.abs(lu[pivot][k])<=Number.EPSILON*scale*n)throw Error('Numerically singular LU pivot at '+k);[lu[k],lu[pivot]]=[lu[pivot],lu[k]];[permutation[k],permutation[pivot]]=[permutation[pivot],permutation[k]];for(let i=k+1;i<n;i++){lu[i][k]/=lu[k][k];for(let j=k+1;j<n;j++)lu[i][j]-=lu[i][k]*lu[k][j];}if(trace)steps.push({k,pivot,lu:lu.map(r=>r.slice()),permutation:permutation.slice()});}return {lu,permutation,steps};}
function solve(F,b){const {lu,permutation}=F,n=b.length,y=permutation.map(i=>b[i]),x=Array(n).fill(0);for(let i=0;i<n;i++)for(let j=0;j<i;j++)y[i]-=lu[i][j]*y[j];for(let i=n-1;i>=0;i--){let r=y[i];for(let j=i+1;j<n;j++)r-=lu[i][j]*x[j];x[i]=r/lu[i][i];}return {x,y,permutedRhs:permutation.map(i=>b[i])};}
const transpose=A=>A[0].map((_,j)=>A.map(r=>r[j]));
function adjoint(F,u,trace=false){const A=Array.from({length:F.N},()=>Array(F.N).fill(0));F.matrix(u,false).forEach((row,i)=>row.forEach(([j,v])=>A[i][j]=v));const AT=transpose(A),rhs=F.dJdu(u).map(v=>-v),fac=factor(AT,trace),s=solve(fac,rhs),res=F.vjp(u,s.x).map((v,i)=>v-rhs[i]),max=a=>Math.max(...a.map(Math.abs)),den=Math.max(...AT.map(r=>r.reduce((s,v)=>s+Math.abs(v),0)))*max(s.x)+max(rhs);return {...s,A,AT,rhs,fac,residual:K.norm(res),backwardError:den?max(res)/den:0,gradient:F.gradient(u,s.x)};}
const api={factor,solve,transpose,adjoint};if(typeof module!=='undefined')module.exports=api;else root.BDirect=api;
})(typeof window!=='undefined'?window:null);
附录 D:光滑性、limiter 与求导收敛
Appendix D: smoothness, limiters and derivative convergence
先分辨四种不同问题:原始残差是否收敛、经典导数是否存在、固定线性系统是否解收敛、网格加密后是否逼近物理敏感度。实验不预设硬 limiter 必然发散,也不把光滑 limiter 当作收敛保证。
Separate four questions: convergence of the primal residual, existence of a classical derivative, convergence of a fixed linear solve, and approximation of physical sensitivities under refinement. A hard limiter need not diverge, and a smooth limiter does not guarantee convergence.
本附录独立于正文无黏源项模型,使用自己的参数和单元平均网格。格点页也显示相同的有限体积参考实验;不能把它误认为格点格式的验证。
This appendix has its own controls and cell-average grid, independent of the main inviscid source model. The node page displays the same finite-volume reference experiment; it is not validation of the node scheme.
D.1 光滑的线性分布,也会触及 minmod 切换
D.1 Even a smooth linear profile can reach a minmod switch
μ=0 的左右导数分别为 +1 与 −1。程序相等时选择左分支,只给出该分支的线性化。中心差分得到零不代表导数存在。以下交互同时列出切换条件、局部 tangent 和反向累加权重。
At μ=0 the left and right derivatives are +1 and −1. The code selects the left branch at a tie, defining that branch’s linearization only. A zero central difference does not establish a derivative. Inspect the switch, local tangent and reverse accumulation weights below.
D.2 完整导数与冻结系数究竟差在哪里
D.2 Full derivatives versus frozen coefficients
冻结 limiter 系数删除最后一项;保存 minmod 的分支仍应对分支内部运算完整求导。两者都可能通过各自的转置点积测试,但冻结系数一般不是原始残差的完整导数。为了明确对比,冻结模式保持基态重构值相同,只替换线性化;它没有重新定义一个全局相同的可微残差。
Freezing the coefficient drops the final term. Saving a minmod branch still differentiates all operations within that branch. Both maps may pass their own transpose test, while coefficient freezing generally is not the full residual derivative. The frozen mode preserves the base reconstruction value and replaces only its linearization; it does not define the same globally differentiable residual.
光滑对照采用上述有理函数,对全部实数输入光滑;不再附加硬 min/max 或符号截断。它保持常值状态,但不是本页证明过的 TVD 或正性保持格式。ε 使用无量纲变量差分的尺度,改变 ε 就改变离散模型。Euler 若产生非正密度或压力、或超过固定耗散波速上界,将拒绝该状态而不静默裁剪。
The rational comparison is smooth for all real inputs, with no added hard min/max or sign clipping. It preserves constant states, but is not claimed to be TVD or positivity preserving. Epsilon uses the scale of nondimensional variable differences; changing it changes the discretization. Euler rejects nonpositive states or states beyond its fixed dissipation-speed bound rather than silently clipping them.
逐步推导:minmod 的单侧导数与有理限幅的导数
Step by step: the one-sided derivatives of minmod and the derivatives of the rational limiter
在 D.1 的线性分布上,两侧差为 与 。 时二者同号,minmod 取绝对值较小者:
On the linear profile of D.1 the two differences are and . For they have the same sign and minmod takes the smaller magnitude:
中心差分为零,却既不是左导数也不是右导数;程序在相等时选左分支 ,返回的线性化是 。有理对照把两个差 、 放进一个商,求导用商法则:
The central difference is zero, yet it is neither one-sided derivative; the program picks the left branch at the tie and returns the linearisation . The rational comparison puts the two differences , into one quotient, differentiated by the quotient rule:
在相等点 ,两个偏导相同:, 时趋于 ,也就是中心平均的导数。本页默认 、 时为 。沿 D.1 的切片(两差为 、), 在 按对称性为零:有理函数没有尖点。minmod 在同一点的所选分支导数是 ,与 不同:两个模型在光滑单调区域也给出不同的线性化。
At the tie the two partials coincide: , which tends to , the derivative of the central average, as . With this page's default and it equals . Along the slice of D.1 (differences and ), vanishes at by symmetry: the rational function has no kink. The selected-branch derivative of minmod at the same point is , not : the two models linearise differently even in a smooth monotone region.
D.3 无激波稳态 Burgers:分支、矩阵和误差
D.3 Smooth steady Burgers: branches, matrices and errors
两端给定解析 Dirichlet 值,源项积分由解析通量端点差计算。每次参数扰动都同时改变源项和边界,但数值未知量仍由残差求解,未填入解析解。设计变量为振幅 a,目标 J=hΣ(1+xi)ui。一阶与二阶使用相同 Rusanov 波速 2 和中心黏性通量;边界扩散距离为 h/2,首末单元不重构。
Both ends use analytic Dirichlet data; integrated sources are analytic physical-flux differences. Perturbations change sources and boundaries together, but unknowns are solved from residuals, never filled with the analytic profile. The design variable is amplitude a and J=hΣ(1+xi)ui. Both orders use Rusanov speed 2 and centred viscous fluxes; boundary diffusion distances are h/2, with no endpoint-cell reconstruction.
U = sampled reference as an initial guess only
repeat until ||R(U,a)|| is small:
A = exact selected-branch Jacobian
d = LU(A, -R) # small scalar experiment
line-search U + step*d # reject failed convergence
A v = -R_a; A^T psi = -J_U^T
g = J_a + psi^T R_a
compare full-rerun FD, branch changes and LU/GS/GMRES本标量实验用小矩阵 LU 及回溯得到原始状态,以隔离导数线性系统的求解困难。随后在同一个固定状态上比较 GS 缺陷校正、GS 预处理 GMRES 与 LU。精确算子通过计算图 JVP/VJP;一阶矩阵仅构造预处理,伴随使用整个预处理映射的转置。
This scalar experiment uses small-matrix LU and backtracking for the primal, isolating difficulty in derivative linear solves. At the same fixed state, compare GS defect correction, GS-preconditioned GMRES and LU. Exact actions use graph JVP/VJP; a first-order matrix builds the preconditioner, whose entire map is transposed for the adjoint.
GS 发散而 LU 后向残差很小,不能解释成“导数不存在”。反过来,求解器收敛也不能证明 minmod 切换点存在唯一导数。下方显示线性真残差、无穷范数条件数、迭代矩阵范数、原始分支变化和完整重算差分。没有计算精确谱半径:范数大于 1 不是发散证明。
GS divergence with a small LU backward residual does not establish nonexistence of a derivative. Conversely, linear convergence does not establish a unique derivative at a minmod switch. Inspect true residuals, infinity-norm conditioning, iteration-matrix norm, primal branch changes and full-rerun differences. Exact spectral radii are not computed: a norm above one is not proof of divergence.
D.4 固定物理时间:峰值增长还是求导失稳
D.4 Fixed physical time: concentration or derivative instability
独立非定常实验使用初值的精确单元平均,左/右外侧状态为 1/0,T=0.12。选择 ν、初始过渡层宽度 δ、网格及 limiter;只对位置 a=0.4 求导。SSP-RK2 的时间步数由网格、固定波速 1.5 和黏性稳定性上界确定,在位置差分中保持不变。所有 RK 子阶段的完整导数都进入 tangent 与反向时间 adjoint。
The independent transient experiment uses exact initial cell averages, exterior states 1/0 and T=0.12. Choose viscosity, initial width, grid and limiter; differentiate only position a=0.4. SSP-RK2 step count uses the mesh, fixed speed 1.5 and a viscous stability bound, and remains unchanged in position perturbations. Both RK stages contribute fully to the tangent and reverse-time adjoint.
黏性 ν 控制 PDE,ε 控制重构映射,δ 控制初始层宽度,三者不能混为一谈。先固定 ν、δ 加密,再改变 ν 或 δ。图同时显示 max|v|、hΣ|v| 和积分梯度;对于移动间断的无黏极限,tangent 可以集中为分布,峰值增大而积分敏感度仍有限。这里不给一般激波的连续伴随收敛保证。
Viscosity changes the PDE, epsilon changes reconstruction, and delta changes initial thickness. First refine at fixed viscosity and width, then vary one parameter. Compare max|v|, hΣ|v| and integral gradients. A moving discontinuity’s inviscid-limit tangent can concentrate as a distribution while integrated sensitivity stays finite. No general continuous-adjoint convergence guarantee is asserted.
逐步推导:初值导数为什么越来越尖
Step by step: why the initial derivative sharpens
对初值直接求导, 的导数是 :
Differentiate the initial data directly; the derivative of is :
峰值与 成正比,面积却恒为 : 时导数趋于 。离散初值用精确单元平均,单元 的导数是 ,不超过 :峰值由 决定。所以“最大值随加密增大”是分布型导数的正常表现;应当检查的是积分量(例如 )是否收敛。
The peak is proportional to while the area stays : as the derivative tends to . The discrete initial data use exact cell averages, whose derivative in cell is , at most : the peak is set by . A maximum that grows under refinement is therefore normal for a distributional derivative; what must converge are integrated quantities such as .
D.5 阅读结果的顺序
D.5 How to interpret results
先检查原始真残差与状态有效性,再判断 limiter 分支和单侧导数,再比较同一个线性系统的 GS/GMRES/LU,最后进行网格加密及连续参考比较。两种转置矩阵有相同的二范数条件数;迭代差异要从预处理、右端项、非正规增长及停止标准解释,不能笼统归因于“伴随更不稳定”。
Check primal residual and admissibility first, then limiter branches and one-sided derivatives, then GS/GMRES/LU on the same system, and finally refinement against continuous references. A matrix and its transpose have the same 2-norm condition number; iteration differences require examining preconditioners, right-hand sides, nonnormal growth and stopping criteria, not a blanket claim that adjoints are less stable.
查看完整独立标量实验与线性诊断实现Complete scalar experiment and linear diagnostics
/* Small educational systems: matrix-free exact products; explicit matrices only for LU/conditioning diagnostics. */
(function(root){'use strict';
const Q=typeof module!=='undefined'?require('../q1d/reconstruction/math'):root.QMath;
const {norm,dot}=Q,zeros=n=>Array(n).fill(0),transpose=A=>A[0].map((_,j)=>A.map(r=>r[j])),mv=(A,x)=>A.map(r=>dot(r,x));
function factor(A){const n=A.length,lu=A.map(r=>r.slice()),p=Array.from({length:n},(_,i)=>i),scale=Math.max(...A.flat().map(Math.abs));for(let k=0;k<n;k++){let j=k;for(let i=k+1;i<n;i++)if(Math.abs(lu[i][k])>Math.abs(lu[j][k]))j=i;if(!Number.isFinite(lu[j][k])||Math.abs(lu[j][k])<1e-14*scale)throw Error('Singular LU pivot '+k);[lu[j],lu[k]]=[lu[k],lu[j]];[p[j],p[k]]=[p[k],p[j]];for(let i=k+1;i<n;i++){lu[i][k]/=lu[k][k];for(let j=k+1;j<n;j++)lu[i][j]-=lu[i][k]*lu[k][j];}}return {lu,p};}
function solve(F,b){const {lu,p}=F,n=b.length,x=p.map(i=>b[i]);for(let i=0;i<n;i++)for(let j=0;j<i;j++)x[i]-=lu[i][j]*x[j];for(let i=n-1;i>=0;i--){for(let j=i+1;j<n;j++)x[i]-=lu[i][j]*x[j];x[i]/=lu[i][i];}return x;}
function assemble(L,n){const A=Array.from({length:n},()=>zeros(n));for(let j=0;j<n;j++){const e=zeros(n);e[j]=1;L.jvp(e).R.forEach((v,i)=>A[i][j]=v);}return A;}
function limiter(a,b,kind='minmod',eps=.01){if(kind==='first')return {s:0,da:0,db:0,branch:'first',margin:1};if(kind==='none')return {s:(a+b)/2,da:.5,db:.5,branch:'unlimited',margin:1};if(kind==='smooth'){const den=a*a+b*b+eps*eps,num=a*b*(a+b);return {s:num/den,da:((2*a*b+b*b)*den-2*a*num)/(den*den),db:((a*a+2*a*b)*den-2*b*num)/(den*den),branch:'smooth',margin:1};}const margin=Math.min(Math.abs(a),Math.abs(b),Math.abs(Math.abs(a)-Math.abs(b)));if(a*b<=0)return {s:0,da:0,db:0,branch:'zero',margin};return Math.abs(a)<=Math.abs(b)?{s:a,da:1,db:0,branch:'left',margin}:{s:b,da:0,db:1,branch:'right',margin};}
function slope(o,a,b,kind,eps,mode='full'){const r=limiter(o.value(a),o.value(b),kind,eps),{c,add,mul,sub,div}=o;if(mode==='frozen'){const d=(o.value(a)+o.value(b))/2,phi=Math.abs(d)>1e-15?r.s/d:0;return add(c(r.s),mul(c(phi/2),add(sub(a,c(o.value(a))),sub(b,c(o.value(b))))));}if(kind==='smooth'){const den=add(add(mul(a,a),mul(b,b)),c(eps*eps));return div(mul(mul(a,b),add(a,b)),den);}return add(c(r.s),add(mul(c(r.da),sub(a,c(o.value(a)))),mul(c(r.db),sub(b,c(o.value(b))))));}
function operators(U,p,evalFn,mode='full'){const t=Q.tape();t.op.value=i=>t.val[i];const ui=U.map(t.op.input),pi=p.map(t.op.input),r=evalFn(t.op,ui,pi,mode);return {R:r.R.map(i=>t.val[i]),J:t.val[r.J],jvp:(v,dp=zeros(p.length))=>{const z=t.forward([...v,...dp]);return {R:r.R.map(i=>z[i]),J:z[r.J]};},vjp:(w,j=0)=>{const z=t.reverse([...r.R,r.J],[...w,j]);return {U:z.slice(0,U.length),p:z.slice(U.length)};}};}
const real={...Q.real,value:x=>x};
function gsMatrix(P,sweeps=2){const n=P.length,B=Array.from({length:n},()=>zeros(n));for(let k=0;k<n;k++){const b=zeros(n),x=zeros(n);b[k]=1;for(let s=0;s<sweeps;s++)for(let i=0;i<n;i++){let v=b[i];for(let j=0;j<n;j++)if(i!==j)v-=P[i][j]*x[j];if(Math.abs(P[i][i])<1e-14)throw Error('Zero GS diagonal');x[i]=v/P[i][i];}x.forEach((v,i)=>B[i][k]=v);}return B;}
function compare(L,P,rhs,adj=false,omega=.7,preAction=null){const n=rhs.length,A=assemble(L,n),AA=adj?transpose(A):A,F=factor(AA),exact=solve(F,rhs),B0=preAction?transpose(Array.from({length:n},(_,j)=>{const e=zeros(n);e[j]=1;return preAction(e);})):gsMatrix(P),B=adj?transpose(B0):B0,op=v=>adj?L.vjp(v).U:L.jvp(v).R,pre=v=>mv(B,v),gm=Q.gmres(op,rhs,pre,1e-11,n);let x=zeros(n),history=[],reason='limit';for(let k=0;k<250;k++){const ax=op(x),r=rhs.map((b,i)=>b-ax[i]),rn=norm(r);history.push(rn);if(rn<1e-10){reason='converged';break;}if(!Number.isFinite(rn)||rn>Math.max(1,norm(rhs))*1e8){reason='growth';break;}const d=pre(r);x=x.map((v,i)=>v+omega*d[i]);}
const inv=transpose(Array.from({length:n},(_,j)=>{const e=zeros(n);e[j]=1;return solve(F,e);})),ni=M=>Math.max(...M.map(r=>r.reduce((s,v)=>s+Math.abs(v),0))),cond=ni(AA)*ni(inv),M=AA.map((_,i)=>AA[0].map((_,j)=>(i===j?1:0)-omega*B[i].reduce((s,v,k)=>s+v*AA[k][j],0)));return {exact,luResidual:norm(op(exact).map((v,i)=>v-rhs[i])),gm,gs:{x,history,reason},conditionInfinity:cond,iterationNormInfinity:ni(M),A};}
const api={...Q,zeros,transpose,mv,factor,solve,assemble,limiter,slope,operators,real,gsMatrix,compare};if(typeof module!=='undefined')module.exports=api;else root.EduNumerics=api;
})(typeof window!=='undefined'?window:globalThis);
/* Independent scalar smoothness experiment. Cell averages, fixed mesh and fixed physical time grid. */
(function(root){'use strict';const N=typeof module!=='undefined'?require('./numerics'):root.EduNumerics,{real,zeros,norm}=N;
function model({n=16,kind='minmod',eps=.01,nu=.01,amplitude=.2,phase=.07,omega=.7}={}){
if(!Number.isInteger(n)||n<8||n>128||nu<0||nu>.1||eps<=0)throw Error('Invalid scalar experiment');const h=1/n,x=Array.from({length:n},(_,i)=>(i+.5)*h),params=[amplitude,phase],alpha=2;
function exact(o,x,p){return o.add(o.c(1),o.mul(p[0],o.sin(o.mul(o.c(2*Math.PI),o.sub(o.c(x),p[1])))));}
function derivative(o,x,p){return o.mul(o.mul(o.c(2*Math.PI),p[0]),o.cos(o.mul(o.c(2*Math.PI),o.sub(o.c(x),p[1]))));}
function evaluate(o,U,p,mode='full',lim=kind){const {c,add,sub,mul}=o,slopes=U.map((u,i)=>i===0||i===n-1?c(0):N.slope(o,sub(u,U[i-1]),sub(U[i+1],u),lim,eps,mode)),F=[];
for(let f=0;f<=n;f++){const L=f?add(U[f-1],mul(c(.5),slopes[f-1])):exact(o,0,p),R=f<n?sub(U[f],mul(c(.5),slopes[f])):exact(o,1,p);const left=f?U[f-1]:L,right=f<n?U[f]:R,distance=f===0||f===n?h/2:h;F[f]=sub(sub(mul(c(.25),add(mul(L,L),mul(R,R))),mul(c(alpha/2),sub(R,L))),mul(c(nu/distance),sub(right,left)));}
const physicalFlux=z=>{const u=exact(o,z,p);return sub(mul(c(.5),mul(u,u)),mul(c(nu),derivative(o,z,p)));};const R=U.map((u,i)=>sub(sub(F[i+1],F[i]),sub(physicalFlux((i+1)*h),physicalFlux(i*h)))),J=U.reduce((s,u,i)=>add(s,mul(c(h*(1+x[i])),u)),c(0));return {R,J};}
function linearize(U,p=params,mode='full',lim=kind){return N.operators(U,p,(o,u,p,m)=>evaluate(o,u,p,m,lim),mode);}
function branches(U){return U.map((u,i)=>i===0||i===n-1?'end':N.limiter(u-U[i-1],U[i+1]-u,kind,eps).branch);}
function solve(p=params,start){let U=start?start.slice():x.map(xx=>exact(real,xx,p)),history=[],switches=[],prev=branches(U),states=[];for(let k=0;k<70;k++){const L=linearize(U,p),r=norm(L.R);history.push(r);states.push(U.slice());const branch=branches(U);switches.push(branch.filter((v,i)=>v!==prev[i]).length);prev=branch;if(r<1e-11)return {U,history,switches,states,converged:true};const d=N.solve(N.factor(N.assemble(L,n)),L.R.map(v=>-v));let t=1,accepted=false;for(let j=0;j<25;j++){const v=U.map((u,i)=>u+t*d[i]);if(Math.min(...v)>.1&&Math.max(...v)<1.9&&norm(evaluate(real,v,p).R)<(1-1e-4*t)*r){U=v;accepted=true;break;}t*=.5;}if(!accepted)return {U,history,switches,states,converged:false,reason:'No residual-decreasing branch correction'};}return {U,history,switches,states,converged:false,reason:'Iteration limit'};}
function analyse(s,mode='full',parameter=0){if(!s.converged)throw Error('Steady primal has not converged');const L=linearize(s.U,params,mode),P=N.assemble(linearize(s.U,params,'full','first'),n),dp=[0,0];dp[parameter]=1;const d=L.jvp(zeros(n),dp),j=L.vjp(zeros(n),1),t=N.compare(L,P,d.R.map(v=>-v),false,omega),a=N.compare(L,P,j.U.map(v=>-v),true,omega),gradient=L.vjp(a.exact,1).p[parameter],tg=L.jvp(t.exact,dp).J,w=x.map(xx=>Math.cos(3*xx)),v=x.map(xx=>Math.sin(7*xx));return {L,t,a,gradient,tangentGradient:tg,dotError:Math.abs(N.dot(w,L.jvp(v).R)-N.dot(v,L.vjp(w).U))};}
return {n,h,x,params,kind,eps,nu,solve,analyse,linearize,branches,evaluate:(U,p=params)=>evaluate(real,U,p),exact:xx=>exact(real,xx,params)};
}
function transient({n=40,nu=.005,eps=.01,kind='minmod',width=.05,T=.12,position=.4}={}){
if(n<8||n>160||nu<0||width<=0||eps<=0||T<0||T>.2)throw Error('Invalid transient experiment');const h=1/n,x=Array.from({length:n},(_,i)=>(i+.5)*h),alpha=1.5,steps=Math.max(1,Math.ceil(T/(.2/(alpha/h+2*nu/(h*h))))),dt=T/steps;
function spatial(o,U,p,mode){const {c,add,sub,mul}=o,sl=U.map((u,i)=>i===0||i===n-1?c(0):N.slope(o,sub(u,U[i-1]),sub(U[i+1],u),kind,eps,mode)),F=[];for(let f=0;f<=n;f++){const L=f?add(U[f-1],mul(c(.5),sl[f-1])):c(1),R=f<n?sub(U[f],mul(c(.5),sl[f])):c(0),ul=f?U[f-1]:c(1),ur=f<n?U[f]:c(0);F[f]=sub(sub(mul(c(.25),add(mul(L,L),mul(R,R))),mul(c(alpha/2),sub(R,L))),mul(c(nu/(f===0||f===n?h/2:h)),sub(ur,ul)));}return {R:U.map((_,i)=>mul(c(1/h),sub(F[i],F[i+1]))),J:c(0),boundaryFlux:sub(F[0],F[n])};}
// Exact cell averages of a tanh transition; initial derivative follows its moving centre.
const logcosh=z=>Math.abs(z)+Math.log1p(Math.exp(-2*Math.abs(z)))-Math.log(2),initial=a=>x.map((_,i)=>.5-.5*width/h*(logcosh(((i+1)*h-a)/width)-logcosh((i*h-a)/width))),initialDerivative=a=>x.map((_,i)=>.5/h*(Math.tanh(((i+1)*h-a)/width)-Math.tanh((i*h-a)/width)));
function run(a=position,keep=true,mode='full'){let U=initial(a),v=initialDerivative(a),history=[U.slice()],vs=[v.slice()],tape=[],massError=0;
for(let k=0;k<steps;k++){const A=N.operators(U,[],spatial,mode),y=U.map((u,i)=>u+dt*A.R[i]),B=N.operators(y,[],spatial,mode),next=U.map((u,i)=>.5*(u+y[i]+dt*B.R[i]));const va=A.jvp(v).R,vy=v.map((z,i)=>z+dt*va[i]),vb=B.jvp(vy).R;v=v.map((z,i)=>.5*(z+vy[i]+dt*vb[i]));massError=Math.max(massError,Math.abs(h*next.reduce((s,z,i)=>s+z-U[i],0)-.5*dt*(spatial(real,U,[],mode).boundaryFlux+spatial(real,y,[],mode).boundaryFlux)));U=next;if(keep){tape.push({A,B});history.push(U.slice());vs.push(v.slice());}}
const weights=x.map(xx=>h*(1+xx)),J=N.dot(weights,U);let w=weights.slice(),dual=Array(steps+1);dual[steps]=w.slice();if(keep)for(let k=steps-1;k>=0;k--){const {A,B}=tape[k],bw=B.vjp(w).U,wy=w.map((z,i)=>.5*(z+dt*bw[i])),aw=A.vjp(wy).U;w=w.map((z,i)=>.5*z+wy[i]+dt*aw[i]);dual[k]=w.slice();}return {U,J,history,vs,dual,forward:N.dot(weights,v),adjoint:keep?N.dot(w,initialDerivative(a)):null,maxTangent:Math.max(...v.map(Math.abs)),l1Tangent:h*v.reduce((s,z)=>s+Math.abs(z),0),massError};}
return {n,h,x,steps,dt,run,position,nu,width};}
const api={model,transient};if(typeof module!=='undefined')module.exports=api;else root.SmoothBurgers=api;
})(typeof window!=='undefined'?window:globalThis);