简介这份资源面向具有计算物理或工程数值仿真背景的研究生及以上研究人员围绕复现论文《High-order central Hermite WENO schemes:Dimension-by-dimension moment-based reconstructions》展开重点解决一维与二维双曲守恒律的高阶数值求解问题。内容以Python为工具系统实现HWENO空间重构并对比Lax-Wendroff与NCE-RK两种时间离散化方法涵盖网格参数设定、通量计算、边界条件处理、误差计算与收敛性测试最终给出较完整的求解器原型。资源包为单个docx文档约25KB内含详细可运行代码片段与逐段解释便于读者边读边调试。目前已有76人学习。通过这份资料读者可掌握HWENO重构中线性权重、光滑度指标与非线性权重的编码思路理解高阶有限体积法在CFD中的落地方式并获得从一维推广到二维的完整实现参考适合希望深入理解高级WENO格式程序设计的数值分析学习者。1. 从 Euler 方程的激波说起为什么一维和二维 HWENO 值得你花时间做计算流体力学的人迟早会撞上同一个问题用高阶格式算激波波前波后总冒出非物理振荡把网格加密到怀疑人生也压不住。WENO 系列格式就是冲着这个矛盾来的——在光滑区保持高阶精度在间断附近自动降阶、抑制振荡。而 HWENOHermite Weighted Essentially Non-Oscillatory在此基础上多做了一件事把解的导数信息也纳入重构用同样的模板宽度换到更高的精度和更紧致的模板。一维问题里它意味着更少的网格点达到同样的分辨率二维问题里它意味着在激波、涡、边界层共存时格式的鲁棒性明显好于纯 WENO。这篇笔记面向已经写过一维 Euler 求解器、想把手上的 WENO 换成 HWENO 的从业者也面向刚学完数值方法、想找一个能真正跑起来、能看到收敛阶的练手项目的人。Python 在这里不是玩具——NumPy 的向量化足够把一维和二维的算例在笔记本上跑到可接受的时间而代码可读性让你能一行行看清重构、通量、时间推进到底在干什么。下面从格式构造讲到可运行代码再到二维扩展和踩坑记录目标是让你照着敲完能亲眼看到 L1 误差随网格加密按预期阶数下降。2. HWENO 重构到底比 WENO 多了什么从 Hermite 插值到权重设计2.1 从 WENO 的短板说起标准 WENO-JS 格式在单元 $I_i$ 上用 $r$ 个候选模板做加权组合每个模板上做 $k$ 次多项式重构光滑区权重趋近理想权重达到 $2k-1$ 阶。问题出在两点一是模板宽度随精度阶数增长很快五点、七点模板在边界附近处理麻烦二是只用了单元平均值没有利用导数信息信息利用率低。HWENO 的核心思路是除了单元平均值 $\bar{u}_i$还把导数近似 $u_i$或 $u_x, u_y$作为额外已知量参与重构。这样每个小模板上可以用 Hermite 插值构造更高次多项式在模板宽度不变的情况下把精度提上去。常见做法是构造两个候选模板一个偏左、一个偏右各自用 Hermite 插值得到重构值再用 WENO 型非线性权重组合。2.2 Hermite 插值重构的构造以一维标量守恒律 $u_t f(u)x 0$ 为例。设已知单元平均值 $\bar{u}i$ 和单元中心导数 $ui$。在模板 $S_0 {I{i-1}, I_i}$ 上用 $\bar{u}{i-1}, \bar{u}i, u{i-1}, ui$ 四个量可以构造三次 Hermite 插值多项式得到 $I{i1/2}$ 处的左重构值 $u^-{i1/2}$。同理在 $S_1 {I_i, I_{i1}}$ 上构造另一个重构。两个候选值做非线性加权$$ u^-{i1/2} \omega_0 , u^{(0)}{i1/2} \omega_1 , u^{(1)}_{i1/2} $$权重 $\omega_k$ 由光滑指示器 $\beta_k$ 决定形式与 WENO-JS 一致$$ \omega_k \frac{\alpha_k}{\sum_j \alpha_j}, \quad \alpha_k \frac{d_k}{(\varepsilon \beta_k)^2} $$其中 $d_k$ 是理想权重由精度阶条件确定。关键区别在于Hermite 插值的理想权重和光滑指示器都要重新推导不能直接套 WENO 的公式。我一般会先用线性权重$\omega_k d_k$验证光滑区的收敛阶确认重构本身没错再打开非线性权重看激波捕捉效果。2.3 导数怎么来HWENO 的闭环问题HWENO 需要导数 $u_i$但导数本身也是未知量。常见处理有两种一是把导数作为独立演化量对 $u_t f(u)_x 0$ 求导得到 $u_t f(u) ux 0$和原方程耦合求解这就是真正的 HWENO 有限体积/有限差分框架二是用相邻单元平均值做紧致差分近似 $ui \approx (\bar{u}{i1} - \bar{u}{i-1})/(2\Delta x)$实现简单但精度和鲁棒性打折。对想快速上手的人我建议先用第二种近似跑通一维确认重构和通量部分正确再升级到耦合求解。下面代码走的是这条路线。2.4 时间推进与通量空间离散完成后半离散格式写成 $d\bar{u}i/dt -(\hat{f}{i1/2} - \hat{f}{i-1/2})/\Delta x$。通量用 Lax-Friedrichs 分裂$\hat{f}{i1/2} \frac{1}{2}(f(u^-) f(u^) - \alpha(u^ - u^-))$$\alpha \max|f(u)|$。时间推进用三阶 TVD Runge-Kutta# SSP-RK3 时间推进CFL 建议取 0.3~0.5 def rk3_step(u, dt, rhs_func): k1 rhs_func(u) u1 u dt * k1 k2 rhs_func(u1) u2 0.75 * u 0.25 * (u1 dt * k2) k3 rhs_func(u2) u_new (1.0/3.0) * u (2.0/3.0) * (u2 dt * k3) return u_newrhs_func返回 $-\partial_x \hat{f}$。dt由 CFL 条件给出dt CFL * dx / max_wave_speed。一维 Euler 方程的最大波速是 $|u| c$二维还要加上 $y$ 方向贡献取两个方向的最大值。3. 一维 HWENO 的 Python 最小实现从重构到收敛阶验证3.1 网格、初值与主循环骨架先搭一维标量问题的骨架用 Burgers 方程 $u_t (u^2/2)_x 0$ 做验证因为它有解析解可以对比又能产生激波。网格用均匀网格边界用周期或 outflow先跑周期。import numpy as np def init_grid(xmin, xmax, N): dx (xmax - xmin) / N xc xmin (np.arange(N) 0.5) * dx # 单元中心 return xc, dx def init_cond_burgers(xc): # 光滑初值后续会陡化成激波 return np.sin(2.0 * np.pi * xc) 0.5 def max_speed_burgers(u): return np.max(np.abs(u))xc是单元中心坐标dx是网格间距。初值用正弦加常数保证激波在计算域内部形成。max_speed_burgers返回全局最大波速用于 CFL 条件。主循环里每步先算导数近似再做 HWENO 重构然后算通量、更新。3.2 导数近似与 Hermite 重构函数导数用中心差分近似边界用单侧差分。重构函数接收单元平均值数组返回每个界面的左右重构值。def compute_derivative(u, dx): # 周期边界下的中心差分近似 up np.roll(u, -1) um np.roll(u, 1) return (up - um) / (2.0 * dx) def hweno_reconstruct_1d(u, du, dx): N len(u) um np.roll(u, 1) up np.roll(u, -1) dum np.roll(du, 1) dup np.roll(du, -1) # 模板 S0: {i-1, i}Hermite 插值到 i1/2 # 三次 Hermite 在右端点的值系数由单元平均值和导数确定 u0 (um u) / 2.0 dx * (dum - du) / 8.0 # 模板 S1: {i, i1} u1 (u up) / 2.0 dx * (du - dup) / 8.0 # 光滑指示器用相邻单元差分的平方 beta0 (u - um) ** 2 dx ** 2 * (du dum) ** 2 beta1 (up - u) ** 2 dx ** 2 * (dup du) ** 2 eps 1e-6 d0, d1 0.5, 0.5 # 理想权重光滑区等权 a0 d0 / (eps beta0) ** 2 a1 d1 / (eps beta1) ** 2 w0 a0 / (a0 a1) w1 1.0 - w0 u_left w0 * u0 w1 * u1 # i1/2 左重构 u_right np.roll(u_left, 1) # i-1/2 的右重构周期平移 return u_left, u_rightu0和u1的表达式来自三次 Hermite 插值在区间端点的取值公式系数 $1/2$ 和 $1/8$ 是推导结果不要随意改。beta是光滑指示器这里用了简化形式实际文献里会有更精细的构造但用于验证收敛阶足够。eps防止除零取 $10^{-6}$ 是常见做法太小会让权重对微小振荡过于敏感太大会污染光滑区精度。d0, d1在光滑区应该由精度阶条件确定这里取等权是因为两个模板对称对奇数阶精度够用。3.3 通量与时间推进def burgers_flux(u): return 0.5 * u ** 2 def rhs_burgers(u, dx): du compute_derivative(u, dx) uL, uR hweno_reconstruct_1d(u, du, dx) fL burgers_flux(uL) fR burgers_flux(uR) alpha max_speed_burgers(u) # Lax-Friedrichs 数值通量 fhat 0.5 * (fL fR - alpha * (uR - uL)) # 界面通量差分注意周期边界用 roll fhat_right fhat fhat_left np.roll(fhat, 1) return -(fhat_right - fhat_left) / dxfhat是 $i1/2$ 界面的数值通量fhat_left是 $i-1/2$ 界面的用np.roll平移得到。alpha取全局最大波速比局部最大值更耗散但更稳初学先用全局。返回的是 $d\bar{u}_i/dt$。3.4 收敛阶验证用光滑解算 L1 误差激波问题没有解析解验证收敛阶要用光滑解。把初值换成 $u_0(x) \sin(2\pi x)$在激波形成前$t0.1$对比数值解和解析解。Burgers 方程此时可用特征线法求解析解或者干脆用非常密的网格做参考解。def l1_error(u_num, u_ref): return np.mean(np.abs(u_num - u_ref)) # 网格序列 Ns [50, 100, 200, 400] errors [] for N in Ns: xc, dx init_grid(0, 1, N) u init_cond_burgers(xc) t, t_end 0.0, 0.1 while t t_end: dt 0.3 * dx / max_speed_burgers(u) dt min(dt, t_end - t) u rk3_step(u, dt, lambda uu: rhs_burgers(uu, dx)) t dt # 参考解用 N3200 的结果插值这里省略插值细节 errors.append(l1_error(u, u_ref_interp)) # 计算收敛阶 for i in range(1, len(errors)): order np.log(errors[i-1] / errors[i]) / np.log(Ns[i] / Ns[i-1]) print(fN{Ns[i]}, order{order:.2f})如果重构和通量都正确你应该看到阶数在 3 左右受导数近似和光滑指示器简化影响可能略低于理论值。如果阶数只有 1 或 2先检查 Hermite 系数符号再检查通量差分方向。如果阶数震荡多半是eps或beta的问题。提示验证收敛阶时先把非线性权重关掉令w0 d0确认线性格式的阶数再打开。这样能把重构误差和权重误差分开定位。4. 二维 HWENO 的扩展维度分裂与模板设计4.1 二维重构的两种路线二维问题 $u_t f(u)_x g(u)_y 0$ 的 HWENO 扩展有两条路。一是维度分裂先沿 $x$ 方向做一维重构和通量再沿 $y$ 方向做一遍每步都是一维 HWENO。实现简单但分裂误差会限制时间精度且对强各向异性问题不够鲁棒。二是真二维重构在二维模板上做双变量 Hermite 插值直接得到界面上的左右状态。精度和鲁棒性更好但模板设计和权重推导复杂得多。我一般先用维度分裂跑通确认二维框架没问题再换真二维。下面代码走维度分裂路线因为它的每一步都能对应到一维代码调试成本低。4.2 二维数组上的方向重构def hweno_reconstruct_2d(u, dx, dy): # 沿 x 方向 dux (np.roll(u, -1, axis1) - np.roll(u, 1, axis1)) / (2.0 * dx) uLx, uRx hweno_reconstruct_1d_axis(u, dux, dx, axis1) # 沿 y 方向 duy (np.roll(u, -1, axis0) - np.roll(u, 1, axis0)) / (2.0 * dy) uLy, uRy hweno_reconstruct_1d_axis(u, duy, dy, axis0) return uLx, uRx, uLy, uRyhweno_reconstruct_1d_axis是把一维重构函数改成沿指定轴操作内部用np.roll的axis参数。二维 Euler 方程需要分别对密度、动量、能量做重构或者对特征变量做重构。对初学者先对守恒量逐分量重构跑通后再考虑特征分解。4.3 二维 Euler 方程的通量与时间推进二维 Euler 方程守恒量 $U [\rho, \rho u, \rho v, E]^T$通量 $F(U)$ 和 $G(U)$ 分别是 $x$ 和 $y$ 方向的。数值通量用二维 Lax-Friedrichs$$ \hat{F}_{i1/2,j} \frac{1}{2}(F(U^-) F(U^) - \alpha_x (U^ - U^-)) $$$\alpha_x \max(|u| c)$$\alpha_y$ 同理。时间推进仍用 SSP-RK3只是rhs里要同时算 $x$ 和 $y$ 方向的通量差分。def rhs_euler_2d(U, dx, dy, gamma1.4): # U shape: (Nx, Ny, 4) uLx, uRx, uLy, uRy hweno_reconstruct_2d(U, dx, dy) # 计算 F(U) 和 G(U)此处省略具体表达式 Fx euler_flux_x(uLx, gamma) Gx euler_flux_x(uRx, gamma) alpha_x max_wave_speed_x(U, gamma) Fhat 0.5 * (Fx Gx - alpha_x * (uRx - uLx)) # y 方向同理 # 差分得到 dU/dt dUdt -(np.roll(Fhat, -1, axis0) - Fhat) / dx \ -(np.roll(Ghat, -1, axis1) - Ghat) / dy return dUdt二维算例建议从光滑的等熵涡开始它有解析解可以验证收敛阶。然后跑双马赫反射或激波-涡干扰看激波捕捉效果。二维计算量比一维大得多Nx200, Ny200在笔记本上跑几千步大概几分钟用 NumPy 的向量化能接受。注意二维维度分裂时$x$ 和 $y$ 方向的通量差分要分别用对应的dx和dy不要混用。我见过有人把dy写成dx结果收敛阶怎么都对不上查了半天。5. 避坑与排查HWENO 实现里最容易翻车的五个地方5.1 收敛阶只有一阶或二阶现象网格加密L1 误差下降很慢算出来的阶数远低于预期。原因最常见的是 Hermite 插值系数符号错了或者导数近似用了单侧差分但没做边界处理。解决先用线性权重、光滑初值、周期边界跑把重构函数单独拿出来在已知多项式上测试看它能否精确重构。如果重构对了再查通量差分方向。5.2 激波附近出现非物理振荡现象激波前后出现高频波动加密网格也不消失。原因光滑指示器beta的构造不够敏感或者eps太小导致权重在间断附近没有及时偏向光滑模板。解决把eps从 $10^{-6}$ 调到 $10^{-3}$ 试试或者改用更精细的beta构造比如加入二阶差分的平方项。如果还不行检查理想权重d0, d1是否满足精度阶条件。5.3 时间步长稍大就发散现象CFL 取 0.5 以上跑几步就 NaN。原因HWENO 的空间离散比 WENO 更紧致但稳定性区域不一定更大尤其是导数近似用中心差分时有效模板宽度增加CFL 限制更严。解决一维先取 CFL0.3二维取 0.2确认稳定后再逐步调大。如果必须用大 CFL考虑隐式时间推进或加滤波。5.4 二维算例跑得极慢现象NxNy100跑 1000 步要十几分钟。原因Python 循环没向量化或者重构函数里用了逐元素操作。解决所有数组操作都用 NumPy 的切片和roll避免for循环。如果还慢把重构和通量函数用numba.jit加速通常能快 10 倍以上。再不够就上 Cython 或直接写 C 扩展。5.5 边界处理导致精度下降现象内部收敛阶正常但加上边界后整体阶数掉下来。原因HWENO 需要导数信息边界处导数近似只能用单侧差分精度低于内部。解决边界附近降阶处理或者用 ghost 单元加高阶单侧重构。如果做周期问题直接用np.roll绕开边界验证阶段先别碰非周期边界。6. 从能跑到好用HWENO 的进阶调参与验证习惯跑通一维和二维之后下一步是让格式真正好用。我自己的习惯是每次改完重构或通量先跑三个标准算例光滑解收敛阶、Burgers 激波、二维等熵涡。三个都过了再上复杂问题。收敛阶用至少四套网格阶数表打印出来低于理论值 0.3 以上就要查。进阶方向有几个。一是把导数作为独立演化量真正实现 HWENO 的耦合框架精度和鲁棒性都会提升但代码量翻倍。二是换用特征分解重构对 Euler 方程的系统问题更稳尤其是激波和接触间断共存时。三是加限制器或保正处理防止密度或压力出现负值。四是把时间推进换成更高效的 SSP-RK 或隐式格式配合 Jacobian-Free Newton-Krylov 求解器。验证方法上除了收敛阶我还会看数值解的总变差TV是否随时间单调不增这是 TVD 性质的一个粗略检验。另外把数值通量在界面上的值打印出来和精确 Riemann 解对比能快速定位是重构错了还是通量错了。参数方面eps我一般取 $10^{-6}$ 到 $10^{-3}$ 之间具体看问题CFL 一维 0.3~0.5二维 0.2~0.4理想权重d0, d1如果模板对称就取等权不对称要重新推导。这些值不是金科玉律但能让你少走弯路。最后说个血泪教训我最早实现 HWENO 时把 Hermite 插值的系数 $1/8$ 写成了 $1/4$收敛阶怎么都上不去查了两天才发现。后来养成了一个习惯每写一个重构公式先用常数、线性、二次多项式各测一遍看它能不能精确重构。这个习惯帮我省了无数个 debug 的夜晚。希望帮到你。本文还有配套的精品资源点击获取