企业官网建设流程全解析

简介：面向计算流体力学初学者与相关课程学习者的MATLAB实践代码包，围绕一维激波管问题提供可运行的数值求解示例，帮助读者理解可压缩Euler方程、黎曼问题及高分辨率格式的基本实现。压缩包共46个.m文件，整体约22KB，按Lax-Friedrichs、Lax-Wendroff、Roe和ENO等格式分目录组织，涵盖主控脚本、数值通量、边界条件、守恒变量与原始变量转换、精确解、精度阶验证和结果绘图等环节，结构紧凑，便于逐项阅读与对照实验。已有3327人学习，说明其在教学与自学场景中具有一定参考价值。读者可通过运行和修改代码，观察不同格式对激波、接触间断和稀疏波的捕捉差异，并掌握Sod问题设置、通量计算、误差统计、收敛阶检验及曲线可视化的完整流程，适合作为课程作业、算法复现和格式入门的辅助材料。

1. 激波管算例的入口：一维可压缩求解器到底能验证什么

激波管问题是验证一维可压缩求解器最常用的算例，ENO重构、Lax-Wendroff格式和Roe近似黎曼解经常被放在同一套代码里对照。这份资源把三者串成一条可运行链路：给定初始间断，推进若干时间步，输出密度、速度、压力曲线，用来看激波、接触间断和膨胀波是否分开且位置合理。它适合正在学计算流体力学、需要做课程设计或准备二次开发求解器的开发者；如果你只想要一个黑盒软件，它不合适。

拿到包后先别急着改格式，先跑默认参数，确认能出图，再决定是补边界条件还是换通量函数。很多问题不是格式本身，而是初始条件、边界和输出时刻对不上。下面按“目录定位—核心模块—运行验证—踩坑排查—二次开发”的顺序拆。

2. 代码包定位与技术选型：从目录到可执行文件的运行链路

这份资源的核心不是某一个函数，而是一条最小可运行链路：网格和初始场生成、界面左右态重构、数值通量计算、时间推进、边界处理、结果落盘、后处理绘图。常见做法是用 Fortran 90 写主循环，参数放在 namelist 或文本文件里，输出 dat 或 csv，再用 Python 或 gnuplot 画图。选择 Fortran 不是怀旧，而是数组运算、编译优化和数值代码的可读性在中小规模一维求解器里很直接；如果你更熟 C++ 或 Python，照着模块边界替换语言也不难，关键是别把“重构—通量—推进”三段揉成一个函数。

2.1 目录结构与文件职责

拿到压缩包后，我一般先列文件，不急着编译。典型结构如下表，文件名可能不同，但职责基本对应。先找到主程序入口和参数文件，再看通量函数和重构函数是否分开；如果所有逻辑都塞在 main 里，二次开发会很痛苦，但作为学习样板仍然能跑。

目录结构背后是模块边界。主程序只负责“什么时候调用谁”，通量函数只负责“给左右态返回数值通量”，重构函数只负责“给单元值返回界面左右态”。这种拆法有一个好处：换格式时只动一个文件，不用在主循环里翻半天。血泪经验是，很多“结果不对”其实是输出列顺序改了，绘图脚本还按旧列读，看起来像激波位置飞了。

2.2 Roe、ENO 与 Lax-Wendroff 在求解器里的分工

一维 Euler 方程的守恒量通常写成 $U=[\rho,\rho u,E]^T$，通量 $F=[\rho u,\rho u^2+p,u(E+p)]^T$。Roe 格式做的是近似黎曼解：用左右态构造 Roe 平均，算特征值和波强，再组装数值通量。ENO 不是通量格式，它负责在单元值基础上做高阶重构，给通量函数提供更准确的界面左右态；如果重构退化成零阶，Roe 就变成一阶迎风，激波和接触间断会被抹得很厉害。Lax-Wendroff 则是另一条路线，用时间泰勒展开和空间中心差分构造二阶格式，代码短，但在强间断附近容易振荡，常被拿来和 Roe 做对照。

选型上，Roe 适合看激波和接触间断的分离，但必须处理熵修正和正性保护；ENO 适合在光滑区提高精度，但模板选择和边界处理写不好会引入额外振荡；Lax-Wendroff 适合做基线，跑得快、代码少，却不适合直接拿来交强间断算例的最终结果。我一般把三者放在同一套网格和同一套时间推进下比较，先固定 CFL 和末时刻，再换通量/重构，这样差异才可解释。

2.3 编译、参数文件与最小运行闭环

先确认编译器。常见做法是用 gfortran 直接编多个源文件，模块文件先编，主程序后编；如果源文件里有 module，编译顺序不能反。下面这条命令是典型的最小闭环，文件名按实际替换。-O2打开优化，-fcheck=all在调试阶段打开数组越界和未初始化检查，跑通后再关掉。

# 先清理旧模块和可执行文件，避免旧 .mod 干扰 rm -f *.mod *.o solver # 先编译模块文件，再编译主程序；顺序反了会出现“找不到模块” gfortran -O2 -fcheck=all -c roe_flux.f90 gfortran -O2 -fcheck=all -c eno_recon.f90 gfortran -O2 -fcheck=all -c lax_wendroff.f90 gfortran -O2 -fcheck=all -c main.f90 gfortran -O2 -fcheck=all -o solver *.o # 运行；参数文件路径按包里实际位置改 ./solver params.nml

这段命令的逻辑是“模块先编、主程序后编、最后链接”。参数说明：-c只编译不链接；-fcheck=all是调试开关，能抓数组越界和未初始化变量；-O2是常用优化级别，若结果出现 NaN，先降回-O0 -g定位。运行后先看屏幕输出有没有“step、time、dt、min rho、max mach”这类监控量，再看 output 目录有没有新文件。最小验收标准不是“没报错”，而是：时间推进到指定末时刻、密度和压力没有负值、输出文件列数稳定、绘图能看出激波和接触间断。

2.4 输出约定与第一次验收标准

输出文件常见两种：每个时刻一个 dat，或者一个文件里按块分段。无论哪种，先确认列顺序，常见是x rho u p e，也有x rho u p mach。第一次跑不要直接看漂亮图，先做三个检查：第一，初始时刻的密度台阶是否在预期位置；第二，末时刻的激波是否向右、膨胀波是否向左、接触间断是否在两者之间；第三，网格加密一倍后波形是否更锐，而不是出现新振荡。若第一项就不对，问题在初始条件或坐标；若第二项方向反了，检查通量符号和边界；若第三项越加密越差，多半是 CFL 或熵修正。

这个阶段不要急着调格式参数。先把默认算例跑出可重复的结果，记录网格数、CFL、末时刻、输出列和绘图脚本版本。后面换 ENO 阶数或换 Roe 熵修正时，才有对照基线。常见做法是保留一份baseline/输出，不改名、不覆盖，二次开发时每次对比都拿它当参照。

3. 核心模块拆解：Roe通量、ENO重构与时间推进的代码骨架

这一章按数据流拆：先有网格和守恒量，再有界面左右态，再有数值通量，最后时间推进。下面代码是按常见 Fortran 90 写法整理的核心骨架，不是包内原文的逐行复制；读的时候重点看数组下标、左右态顺序和特征值处理，这些地方最容易翻车。

3.1 网格、初始间断与守恒量初始化

一维激波管通常取 $x\in[0,1]$，膜片在 $x=0.5$。Sod 算例左态高密度高压、右态低密度低压，速度初值为零。网格用单元中心存储，界面在单元之间。初始化时把守恒量填进U(1:nx,1:3)，顺序是密度、动量、总能量。总能量由 $E=p/(\gamma-1)+0.5\rho u^2$ 得到。

! 初始化一维激波管：单元中心网格 + Sod 左右态 subroutine init_sod(nx, xmin, xmax, gamma, U, xc, dx) implicit none integer, intent(in) :: nx real(8), intent(in) :: xmin, xmax, gamma real(8), intent(out) :: U(nx,3), xc(nx), dx integer :: i real(8) :: xmid, rhoL, uL, pL, rhoR, uR, pR, EL, ER dx = (xmax - xmin) / real(nx,8) xmid = 0.5d0 * (xmin + xmax) rhoL = 1.0d0; uL = 0.0d0; pL = 1.0d0 rhoR = 0.125d0; uR = 0.0d0; pR = 0.1d0 EL = pL / (gamma - 1.0d0) + 0.5d0 * rhoL * uL * uL ER = pR / (gamma - 1.0d0) + 0.5d0 * rhoR * uR * uR do i = 1, nx xc(i) = xmin + (real(i,8) - 0.5d0) * dx if (xc(i) < xmid) then U(i,1) = rhoL U(i,2) = rhoL * uL U(i,3) = EL else U(i,1) = rhoR U(i,2) = rhoR * uR U(i,3) = ER end if end do end subroutine init_sod

逻辑说明：xc存单元中心，dx是均匀步长，U(i,1:3)分别是密度、动量、总能量。参数说明：nx常用 100～400 做快速验证，gamma空气取 1.4，xmin/xmax取 0/1，膜片位置取中点。注意xmid比较用<还是<=会影响一个单元，但不会改变整体波形；真正要小心的是总能量公式，漏掉动能项会让压力初值直接错。

3.2 Roe平均与数值通量

Roe 通量的关键是 Roe 平均和特征值。给左右态UL、UR，先算 Roe 平均速度、焓、声速，再算三个特征值 $u-a$、$u$、$u+a$，用波强组装通量。强激波附近会出现小特征值，必须加熵修正，否则会出现非物理膨胀激波。下面代码用 Harten 型熵修正，delta是修正阈值，常取 0.1～0.2 倍当地声速。

! Roe 近似黎曼通量：输入左右守恒量，返回界面数值通量 subroutine roe_flux(UL, UR, gamma, delta, F) implicit none real(8), intent(in) :: UL(3), UR(3), gamma, delta real(8), intent(out) :: F(3) real(8) :: rhoL, uL, pL, HL, rhoR, uR, pR, HR real(8) :: sqL, sqR, rhoRoe, uRoe, HRoe, aRoe real(8) :: lam1, lam2, lam3, d1, d2, d3 real(8) :: alpha1, alpha2, alpha3 real(8) :: FL(3), FR(3), drho, du, dp rhoL = UL(1); uL = UL(2)/rhoL pL = (gamma - 1.0d0) * (UL(3) - 0.5d0 * rhoL * uL * uL) HL = (UL(3) + pL) / rhoL rhoR = UR(1); uR = UR(2)/rhoR pR = (gamma - 1.0d0) * (UR(3) - 0.5d0 * rhoR * uR * uR) HR = (UR(3) + pR) / rhoR sqL = sqrt(rhoL); sqR = sqrt(rhoR) rhoRoe = sqL * sqR uRoe = (sqL * uL + sqR * uR) / (sqL + sqR) HRoe = (sqL * HL + sqR * HR) / (sqL + sqR) aRoe = sqrt((gamma - 1.0d0) * (HRoe - 0.5d0 * uRoe * uRoe)) lam1 = uRoe - aRoe; lam2 = uRoe; lam3 = uRoe + aRoe ! Harten 熵修正：避免小特征值导致非物理膨胀激波 lam1 = sign(max(abs(lam1), delta), lam1) lam2 = sign(max(abs(lam2), delta), lam2) lam3 = sign(max(abs(lam3), delta), lam3) drho = rhoR - rhoL; du = uR - uL; dp = pR - pL alpha1 = 0.5d0 * (dp - rhoRoe * aRoe * du) / (aRoe * aRoe) alpha2 = drho - dp / (aRoe * aRoe) alpha3 = 0.5d0 * (dp + rhoRoe * aRoe * du) / (aRoe * aRoe) FL(1) = rhoL * uL FL(2) = rhoL * uL * uL + pL FL(3) = uL * (UL(3) + pL) FR(1) = rhoR * uR FR(2) = rhoR * uR * uR + pR FR(3) = uR * (UR(3) + pR) d1 = abs(lam1) * alpha1 d2 = abs(lam2) * alpha2 d3 = abs(lam3) * alpha3 F(1) = 0.5d0*(FL(1)+FR(1)) - 0.5d0*(d1 + d2 + d3) F(2) = 0.5d0*(FL(2)+FR(2)) - 0.5d0*(d1*(uRoe-aRoe) + d2*uRoe + d3*(uRoe+aRoe)) F(3) = 0.5d0*(FL(3)+FR(3)) - 0.5d0*(d1*(HRoe-uRoe*aRoe) + d2*0.5d0*uRoe*uRoe + d3*(HRoe+uRoe*aRoe)) end subroutine roe_flux

逻辑说明：先由守恒量反算速度和压力，再构造 Roe 平均；三个特征值经过熵修正后参与耗散项。参数说明：delta不能取 0，取 0 会退化成无熵修正的 Roe，强激波附近容易出问题；取太大又会过度耗散，接触间断变糊。我一般先用delta=0.1*aRoe跑 Sod，看压力是否出现负值或膨胀激波，再微调。注意F(3)中中间波的耗散项用0.5*uRoe*uRoe是常见写法，别把总焓项和声速项抄错。

3.3 ENO重构：界面左右态怎么取

ENO 重构的目标是给每个界面提供左右两侧的高阶近似。以密度为例，先看三个候选模板的二阶差分，选绝对值最小的那个，再用对应模板外推。下面是一个简化的三点 ENO 选择骨架，适合二维数组按单元循环调用。实际包里可能按原始变量重构，也可能按特征变量重构；按原始变量实现简单，强激波附近更依赖限制器和熵修正。

! 简化 ENO 重构：给界面 i+1/2 返回左态和右态（以密度为例） subroutine eno_recon_rho(rho, nx, i, rhoL, rhoR) implicit none integer, intent(in) :: nx, i real(8), intent(in) :: rho(nx) real(8), intent(out) :: rhoL, rhoR real(8) :: d1, d2, d3, slope ! 候选模板的差分；越靠边界越要保证下标不越界 d1 = rho(i) - rho(i-1) d2 = rho(i+1) - rho(i) d3 = rho(i+2) - rho(i+1) ! 选绝对值最小的差分作为斜率，等价于选择最光滑模板 if (abs(d1) <= abs(d2) .and. abs(d1) <= abs(d3)) then slope = d1 else if (abs(d2) <= abs(d3)) then slope = d2 else slope = d3 end if rhoL = rho(i) + 0.5d0 * slope rhoR = rho(i+1) - 0.5d0 * slope end subroutine eno_recon_rho

逻辑说明：d1/d2/d3是三个候选模板的斜率，选绝对值最小者可以避开跨间断的模板；rhoL是界面左侧值，rhoR是右侧值。参数说明：i必须从 2 到nx-2循环，边界单元要单独外推或降阶，否则rho(i+2)越界。这个骨架只重构了密度，实际 Euler 系统要对密度、动量、能量分别重构，或者先转原始变量再重构。注意左右态顺序：Roe 通量里UL是左态、UR是右态，反了会出现反向波。

3.4 二阶时间推进与CFL约束

时间推进常用二阶 TVD Runge-Kutta。先算 $L(U)$，再走两步：$U^{(1)}=U^n+\Delta t L(U^n)$，$U^{n+1}=0.5U^n+0.5(U^{(1)}+\Delta t L(U^{(1)}))$。CFL 数控制稳定步长，一维 Euler 的谱半径是 $|u|+a$，取全局最大值后乘 CFL 和 $\Delta x$。下面代码是主循环里的步长和推进骨架，rhs是空间算子，内部按界面循环调用 Roe 通量和 ENO 重构。

! 二阶 RK 推进 + CFL 步长；U 为守恒量，rhs 为空间残差 subroutine advance_rk2(U, nx, dx, gamma, cfl, t_end, t, dt_out, out_step) implicit none integer, intent(in) :: nx, out_step real(8), intent(inout) :: U(nx,3), t real(8), intent(in) :: dx, gamma, cfl, t_end real(8), intent(out) :: dt_out real(8), allocatable :: U1(:,:), L(:,:) real(8) :: dt, smax, a, u, rho, p integer :: i, step allocate(U1(nx,3), L(nx,3)) step = 0 do while (t < t_end) ! 用当前 U 估算最大波速，得到稳定步长 smax = 0.0d0 do i = 1, nx rho = U(i,1) u = U(i,2)/rho p = (gamma-1.0d0)*(U(i,3)-0.5d0*rho*u*u) a = sqrt(max(gamma*p/rho, 1.0d-12)) smax = max(smax, abs(u)+a) end do dt = cfl * dx / max(smax, 1.0d-12) if (t + dt > t_end) dt = t_end - t call spatial_operator(U, L, nx, dx, gamma) U1 = U + dt * L call spatial_operator(U1, L, nx, dx, gamma) U = 0.5d0*U + 0.5d0*(U1 + dt*L) t = t + dt step = step + 1 if (mod(step, out_step) == 0) call write_output(U, nx, t) end do dt_out = dt deallocate(U1, L) end subroutine advance_rk2

逻辑说明：先由当前守恒量算最大波速，得到dt；然后两次调用空间算子完成 RK2。参数说明：cfl常用 0.3～0.5，取 1.0 附近很容易在激波附近炸；out_step控制输出频率，太小会生成大量文件；t_end要和绘图脚本里的时刻一致。注意spatial_operator内部要先做 ENO 重构，再逐界面算 Roe 通量，最后组装残差；如果先算通量再重构，顺序就反了。

4. 边界、输出与验证：把数值结果落到密度曲线上

核心模块跑通后，下一步不是继续加格式，而是让结果可验证。激波管算例的验收点很具体：激波位置、接触间断位置、膨胀波扇面、压力平台。边界条件、输出列和绘图坐标任何一处不一致，都会让人误判格式有问题。我一般先把输出固定成“单元中心 + 密度/速度/压力”的纯文本，再用短脚本画图，图对了再谈加密。

4.1 边界条件与通量循环

一维激波管常用零梯度或固定远场边界。简单实现是把边界外一层单元复制邻近内点，或者直接给边界界面赋左右态。下面代码是零梯度边界的常见写法：在每次空间算子调用前，把U(1,:)和U(nx,:)的外推值填到 ghost 单元。若代码不用 ghost 单元，而是在界面循环里特判，逻辑等价，但更容易把左右态顺序写反。

! 零梯度边界：把内点复制到两端 ghost 单元 subroutine apply_bc_zero_grad(U, nx, ng) implicit none integer, intent(in) :: nx, ng real(8), intent(inout) :: U(1-ng:nx+ng,3) integer :: i, k do k = 1, 3 do i = 1, ng U(1-i,k) = U(1,k) U(nx+i,k) = U(nx,k) end do end do end subroutine apply_bc_zero_grad

逻辑说明：ng是 ghost 层数，ENO 重构用到i+2时至少需要 2 层；边界复制只保证不越界，不保证物理上完全无反射。参数说明：nx是内点数量，数组下标从1-ng到nx+ng；如果主数组只开1:nx，这段要改成界面特判。注意零梯度边界在激波管短时间算例里够用，但末时刻波系接近边界时要加长计算域或换远场边界。

4.2 输出文件格式与Python绘图

输出建议用空白分隔的纯文本，第一行写列名，后面每行一个单元。列顺序固定为x rho u p，不要中途改。下面 Python 脚本读单个时刻文件并画密度曲线，横轴取单元中心。脚本里没有花哨样式，重点是列索引和时刻标注，方便叠图对比。

# 读取一维激波管输出并绘制密度曲线 import numpy as np import matplotlib.pyplot as plt # 文件列顺序：x rho u p；若包内列顺序不同，改 usecols data = np.loadtxt("output/sod_t0.20.dat", comments="#") x = data[:, 0] rho = data[:, 1] p = data[:, 3] plt.figure(figsize=(6, 4)) plt.plot(x, rho, label="rho") plt.plot(x, p, label="p", linestyle="--") plt.xlabel("x") plt.ylabel("value") plt.title("Sod shock tube at t=0.20") plt.legend() plt.tight_layout() plt.savefig("sod_density.png", dpi=150)

逻辑说明：loadtxt读空白分隔文件，comments="#"跳过注释行；usecols没写，默认读全部列，所以列顺序必须和文件一致。参数说明：data[:,1]是密度，data[:,3]是压力；如果输出是x rho u p e，压力仍在第 4 列。注意绘图前先print(data.shape)和print(data[:3,:])，确认没有把时刻列或网格编号列读进来。

4.3 网格加密与结果验收

网格加密不是只看图更细，而是看三个物理特征是否稳定：激波位置随加密变化很小，接触间断不被抹到和激波粘在一起，膨胀波扇面没有出现台阶状振荡。常见做法是跑nx=100、200、400三组，固定 CFL 和末时刻，把密度曲线叠在一张图里。若 100 和 200 差异很大、200 和 400 仍明显分叉，先查边界和 CFL，不要急着把 ENO 阶数调高。

验收时记录四个量：激波前端位置、接触间断位置、膨胀波左端位置、中心压力平台值。位置允许随网格有少量移动，但趋势要一致；压力平台不能出现负值。若曲线在激波后出现高频振荡，多半是通量格式耗散不足或重构跨间断；若接触间断被抹成一个大斜坡，先看 ENO 模板是不是退化成了一阶。

5. 避坑与排查：激波管代码跑不出正确波形时先查这五处

这类代码的报错往往不直接指向根因：编译能过、运行不崩、图也能出，但波形就是不对。下面五条是我在拆一维激波管代码时最常遇到的坑，按“现象—原因—解决”写，先查编译和启动，再查结果形态。

5.1 编译与启动阶段的异常

现象一：编译报“找不到模块”或链接时未定义引用。原因通常不是代码语法，而是编译顺序：use roe_module的主程序先编了，模块文件还没生成；或者旧.mod文件混在目录里，编译器读到了不匹配的接口。解决方法是先删*.mod和*.o，按“模块文件→主程序→链接”的顺序重编；如果模块之间有依赖，把use关系画出来，从最底层模块开始编。不要用通配符一把编所有文件，顺序不可控。

现象二：程序启动即报数组越界或NaN，但 CFL 看起来不大。原因常见于 ghost 层不够：ENO 重构访问i+2，边界只开了 1 层；或者初始场在边界处没填满，U(0,:)是未初始化值。解决方法是把数组下标范围打印出来，确认1-ng到nx+ng都赋了值；调试阶段加-fcheck=all -ffpe-trap=invalid,zero,overflow，让浮点异常直接停下。注意-ffpe-trap在不同编译器上选项名可能不同，先用-fcheck=all兜底。

现象三：输出文件生成了，但绘图脚本读出来只有一列或列数错位。原因通常是输出格式从空格分隔改成了逗号分隔，或者第一行列名被loadtxt当成数据；也可能是不同时刻文件列顺序不一致。解决方法是固定输出格式，第一行写# x rho u p，所有时刻用同一个写文件子程序；绘图前打印data.shape和前三行。不要在主循环里临时改write格式，后面叠图时会对不上。

5.2 结果形态与物理量的异常

现象四：激波位置明显偏左或偏右，但波形形状还行。原因可能是时间步累计误差、输出时刻不对，或者边界条件把波反射回来了。先检查t_end和输出文件名里的时刻是否一致；再检查 CFL 是否过大导致波速被低估；最后看计算域是否太短，激波已经到边界。解决方法是把域长加倍跑一次，若位置回到预期，就是边界反射；若仍偏，检查通量符号和界面左右态顺序。Roe 通量里UL/UR反了，波会朝反方向走。

现象五：接触间断被抹得几乎看不见，或者激波后出现非物理振荡。接触间断抹平常见于 ENO 模板选择过耗散、重构按守恒量直接做、或 Roe 熵修正delta取太大；激波后振荡则常见于熵修正不足、CFL 接近上限、或者 Lax-Wendroff 对照格式被误当成最终格式。解决方法是先把delta调到0.1*aRoe附近跑对照，再把 ENO 重构改成按原始变量做，并检查模板选择是否总选到跨间断的差分。若仍然振荡，降 CFL 到 0.3 跑一次；若振荡消失，问题在时间推进稳定性，不在通量本身。

排查时建议按固定顺序：先看初始时刻输出，再看中间时刻，最后看末时刻。初始时刻不对，查初始化和坐标；中间时刻出现负压，查 CFL 和熵修正；末时刻波系位置不对，查边界和输出时刻。每次只改一个变量，改完保留输出，别在一次运行里同时换网格、换格式、换边界，否则出了问题无法归因。

6. 二次开发：换格式、加限制器与自动化回归的落地做法

这份代码包最适合拿来做“可对照”的二次开发，而不是一上来堆高阶格式。改造顺序我一般这样排：先固定 Sod 算例和输出格式，再把通量函数抽成可替换接口，然后加一个 Lax-Wendroff 对照分支，最后才考虑 ENO 阶数、限制器和特征变量重构。每次替换只动一个模块，主循环和输出不动，这样波形差异才能归因到格式本身。

验证脚本不用复杂，读两个时刻的密度曲线，算激波位置和接触间断位置的相对差，再和基线比。下面这段只做一件事：把baseline和new两条密度曲线插值到同一网格，算最大绝对差，超过阈值就打印警告。它不能替代物理判断，但能防止“图看起来差不多”导致的误判。

# 对比基线和新结果的密度曲线，给出最大绝对差 import numpy as np base = np.loadtxt("baseline/sod_t0.20.dat") new = np.loadtxt("new/sod_t0.20.dat") x = base[:, 0] rho_b = base[:, 1] rho_n = np.interp(x, new[:, 0], new[:, 1]) err = np.max(np.abs(rho_n - rho_b)) print("max |rho_new - rho_base| =", err) if err > 0.05: print("警告：差异超过阈值，先查通量/重构/边界，不要继续加格式")

参数说明：np.interp要求新结果覆盖基线 x 范围；阈值 0.05 只是 Sod 算例的粗筛，网格和末时刻变了要重设。注意别拿不同末时刻的文件对比，文件名里的时刻必须一致。从那以后我每次改通量函数都先把 Sod 问题跑一遍，把密度曲线和基线叠图，再动别的；没有这张图，后面加限制器就是盲调。希望帮到你。

本文还有配套的精品资源，点击获取