简介基于Fortran的二维浅水方程求解器typhon-solver-0.3.0完整源码包面向流体力学数值模拟学习者和研究者可用于复杂地形下的洪水模拟、海洋动力学等场景。软件采用非结构网格上的有限体积法较规则网格能更好适应不规则边界与地形变化。源码包共378个文件以318个f90源文件为主辅以32个make构建脚本及配置文件、说明文档、版本变更记录等整体仅766KB结构精炼便于通读。已有504人学习下载。通过学习源码可系统掌握非结构三角形网格生成、高阶数值积分、Runge-Kutta时间步进、广义拉格朗日乘子法GCL及自由/固定/滑移边界条件处理等核心算法并借助内置山洪暴发、河流流动等算例快速上手适合想深入理解二维浅水方程数值求解或参考Fortran实现高效非结构网格算法的开发者包内还涉及网格提取、插值与CGNS接口等扩展模块便于按需二次开发。1. typhon-solver 源码包在算什么二维浅水非结构网格一个文件名里同时带上 Fortran、二维浅水方程、非结构网格其实已经圈定了用户画像需要模拟河道、潮滩或城市洪涝又不想长期停在“调商业软件”层面的工程师和研究人员。typhon-solver-0.3.0-sources.tar.gz这类源码发布包压缩的是一整套离散算法、网格约定和参数方案而不是几个能直接双击的运行程序。本文要顺着这个源码包把链路拆开先说明非结构网格上浅水方程是怎么离散的再讲怎样把 Fortran 源码编译成可执行文件随后用一个典型算例把边界条件和稳定性参数讲清楚最后落到如何把它接入自己的数据流。适合已经看得懂一两个 Fortran 计算程序、想在浅水模型方向快速建立工程感的读者。2. 非结构网格上的二维浅水方程离散2.1 守恒变量、控制体与数值通量二维浅水方程是沿水深积分后的浅水方程最常见的守恒型写法保留三个守恒变量水深h、x 方向单宽流量qx h*u、y 方向单宽流量qy h*v。写成向量形式为∂U/∂t ∇·F(U) S其中U [h, qx, qy]^TF是包含对流项和静水压力项的通量张量S是源项典型成分包括地形坡度引起的重力投影项和底摩擦项。非结构网格上最匹配的做法是有限体积法把计算域剖成三角形或四边形控制体在每个控制体上做体积分再用高斯定理把面积分转化为边界通量积分。typhon-solver 这类 Fortran 程序的性能差异往往不在方程形式而在网格遍历和通量计算的缓存友好程度。下面这段代码是控制体更新的骨架核心不是复杂而是把“遍历边、取相邻单元、累加通量、除以面积”四步写清楚subroutine update_cell(cell_id, dt, U, F, nfaces, face_list) use grid_data, only: edge_normal, edge_length, cell_area implicit none integer, intent(in) :: cell_id real, intent(in) :: dt real, intent(inout) :: U(:, :) real, intent(in) :: F(:, :, :) ! F(4, nfaces, 2), 通量分量 integer, intent(in) :: nfaces(:), face_list(:, :) integer :: iface, fid real :: flux_sum(4) flux_sum 0.0 do iface 1, nfaces(cell_id) fid face_list(iface, cell_id) flux_sum flux_sum F(:, fid, 1) * edge_normal(1, fid) F(:, fid, 2) * edge_normal(2, fid) end do U(:, cell_id) U(:, cell_id) - dt * flux_sum / cell_area(cell_id) end subroutine update_cell这里把通量F沿外法线方向投影后再累加dt / cell_area放在最后乘是为了避免每个边重复做除法。非结构网格里一个单元少则三条边、多则七八条边边表遍历的顺序会明显影响内存访问效率正式实现里通常会把face_list按单元重排成连续内存块也就是所谓“压缩稀疏行”结构。出于工程判断我在接触新的 Fortran 浅水求解器时会先确认它在边表上是否用了单向遍历而不是每条边重复处理两次。若每条边在两个相邻单元里都被算一遍CPU 开销会翻倍并且容易出现通量不一致导致的局部质量不守恒。2.2 边界通量的黎曼近似HLL、HLLC 怎么选有限体积方法绕不开“边界通量怎么算”的问题。相邻两个单元在同一界面两侧会有不同的守恒变量值界面上的通量要通过近似黎曼解得到。二维浅水模型中最常见的三个方案是 HLL、HLLC 和 Roe。三者的差别用一张表可以看得很清楚格式能否捕捉接触间断干湿边界表现相对计算量典型使用场景HLL不能一般低教学算例、粗网格摸底Roe能容易产生负水深中缓变流水位模拟HLLC能稳定能维持干湿过渡中洪水波、溃坝、复杂地形HLLC 在 HLL 基础上增加了一个接触波估计在浅水方程里表现为对“两个单元间速度差异”的修正。实现时以特征波速做分段多项式重构! HLLC 界面通量的核心选择逻辑 if (sL 0.0) then F_face F_L else if (sL 0.0 .and. sM 0.0) then F_face F_L sL * (U_Lstar - U_L) else if (sM 0.0 .and. sR 0.0) then F_face F_R sR * (U_Rstar - U_R) else F_face F_R end if这段选择逻辑的关键在波速sL、sM、sR。sM是接触波速度sL和sR分别对应左右两侧的波前速度。对于亚临界流动sL为负、sR为正界面通量由中间两个分支决定对于超临界流动两个波速同号通量完全取上游侧。实际编码时要注意干床单元的波速修正否则水深为零的单元会出现除零。我一般会优先保留 HLLC 作为默认通量格式因为它在溃坝、潮滩淹没这类强动边界问题里比 Roe 更皮实。若源码包里同时提供两种格式先用 HLLC 跑通基准算例再用 Roe 做对比算是比较稳妥的验证思路。2.3 时间推进与 CFL 控制显式时间推进在非结构网格浅水模型里用得最多因为它实现简单、天然适合 GPU 并行。稳定性由 CFL 条件约束dt CFL * min_cell(A_cell / Σ λ_edge * L_edge)其中A_cell是控制体面积λ_edge是边上的最大特征波速L_edge是边长。二维非结构网格上特征波速通常取|u| sqrt(g h)其中sqrt(g h)是重力波速。实际 Fortran 代码里这一步往往被包装成一个全局归约dt_local cfl * cell_area(cell_id) / sum_wave_speed(cell_id) call mpi_allreduce(dt_local, dt_global, 1, MPI_REAL, MPI_MIN, comm)对于单 CPU 测试去掉mpi_allreduce直接取最小值即可。CFL 取值在三角形网格上一般取 0.4 到 0.6四边形网格可以放宽到 0.8 左右。低于 0.2 说明网格畸形或者波速计算有误高于 1.0 则基本必发散。3. 从 sources.tar.gz 编译 Fortran 程序的正确姿势3.1 解压、识别构建系统先动手解压源码包用伤害最小的方式判断项目用的什么构建系统tar -xzf typhon-solver-0.3.0-sources.tar.gz cd typhon-solver-0.3.0-sources ls -la file Makefile 2/dev/null ls CMakeLists.txt 2/dev/null解压后重点看三个信号有没有Makefile有没有CMakeLists.txt有没有configure脚本。旧一点的项目经常只有手写 Makefile新一些的会用 CMake。若同时存在优先用 CMake因为它对路径、编译器和外部依赖的管理更干净。源码包内如果出现src/、include/、util/三个目录说明作者已经按模块划分维护起来比单文件大程序容易得多。3.2 把 Fortran 环境装到“能编译”的状态无论目标是调试还是性能测试第一步都是本机的 Fortran 编译器。不同系统安装命令不同但原则一致gfortran配make并行项目再加mpif90。ubuntu 系和 debian 系sudo apt update sudo apt install gfortran gcc make cmake sudo apt install libopenmpi-dev openmpi-bin # 并行版本需要rhel/centos/fedora 系sudo dnf install gcc-gfortran make cmake sudo dnf install openmpi openmpi-devel装完后验证gfortran --version mpif90 --versiongfortran版本差异经常导致ISO_C_BINDING或 OpenMP 语法兼容问题。我通常会专门用一个孤立目录编译源码避免系统里多个编译器版本串扰。若只想快速确认源码能过编译gfortran -c到单个源文件即可无需先链接全部依赖。3.3 Makefile 定制与三个必调编译选项绝大多数 Fortran 数值项目会在 Makefile 顶部分组定义编译器和编译选项。典型的张这样FC gfortran MPIFC mpif90 FFLAGS -O3 -marchnative -cpp -ffree-line-length-none -Wall DBGFLAGS -O0 -g -fcheckall -fbacktrace LIBS -llapack -lblas OBJS main.o swe_mod.o grid_mod.o output_mod.o all: typhon-solver typhon-solver: $(OBJS) $(MPIFC) -o $ $(OBJS) $(LIBS) %.o: %.f90 $(MPIFC) $(FFLAGS) -c $ -o $ debug: FFLAGS $(DBGFLAGS) debug: typhon-solver clean: rm -f *.o *.mod typhon-solver如果是纯串行编译把MPIFC改成gfortran即可。这里有几个选项值得逐一说清楚选项作用什么时候用-O3开启高级自动向量化与内联稳定后的性能测试、生产算例-marchnative按当前 CPU 指令集优化只在本机跑时使用不用于发布-ffree-line-length-none不限制源码行宽处理长行注释和长函数签名-fcheckall运行时检查数组越界、未初始化变量调试阶段必开性能会明显下降-fbacktrace崩溃时打印调用栈定位段错误和浮点异常一个很常见的坑是修改 Makefile 后make不重新编译。因为.mod文件依赖于前一文件但 Makefile 里没有声明模块依赖。若改动grid_mod.f90后其他文件不重编先make clean再make不要手动删.o文件。4. 跑通 typhon-solver 首个算例边界、源项与稳定性参数4.1 用一个矩形水槽做全线摸底第一次拿到可执行文件我不会直接上真实地形而是先用一个规则矩形水槽把“读网格、初始化、推进、输出”四步全部走通。矩形区域 5000 m × 1000 m剖成约 800 个三角形初始水深 1.0 m入口给恒定流量出口固定水位。这样做的理由是解析解容易估算网格质量高发散时问题一定在代码或参数而非网格畸形。控制文件的最简形态我一般写成这样[mesh] file mesh/rect5k_msh.dat [time] t_end 3600 dt_max 0.5 cfl 0.4 [boundary] inlet flux, q2.5 outlet stage, h0.8 wall slip [source] friction manning, n0.025 topography linear, slope0.0001其中inlet给单宽流量q2.5 m²/soutlet给固定水位h0.8 mwall用滑移边界。跑完后先看两个量总水量是否随时间变化以及出入口流量是否达到稳态平衡。若总水量漂移超过 1%多半是通量没有严格守恒问题通常在边界通量实现或地形源项处理。矩形水槽还有一个额外好处可以快速测试网格加密的收敛性。把三角形边长从 100 m 细化到 50 m再加密到 25 m比较同一位置的水位变化。水位差递减且幅度接近 2 倍时说明数值格式达到了二阶精度若水位差没有规律优先检查斜率限制器和干湿阈值。4.2 边界条件与源项的正确打开方式浅水模型的边界条件比普通 CFD 更讲究因为特征波数比 NS 方程少给定变量的个数要严格匹配特征线方向。边界类型可给变量适用场景常见错误入流亚临界流量或水深河道上游、泵站入流同时给水深和流量入流超临界水深、流量急流陡坡只给流量导致特征不完整出流亚临界水位河道下游、海边界给流量导致反射波固壁切向速度堤防、岸线给了法向速度干湿边界无潮间带、洪水漫滩强迫水深大于零地形源项的处理往往比边界条件更影响稳定性。浅水方程中的地形坡度项若采用中心差分静水条件下会因通量梯度和地形源项不平衡而产生虚假流速。常见做法是从通量计算里减去静水压力项或在每个界面上做“水位重构”保证水位 水深 地形高程在静水时处处相等。4.3 三个最容易发散的参数根据我的经验源码包能编译通过却跑不出稳定解的八成都不是编译器问题而是下面三个参数没设好。第一个是干水深阈值eps_dry。低于这个值的水深会被视为干单元不参与通量计算。阈值太大会抹掉真实的浅水流动细节阈值太小会出现水深接近零时的负水深。一般取 0.01 m 到 0.001 m 之间取决于问题尺度。城市洪涝这类小水深大流速的问题阈值要谨慎下调到 0.0001 m。第二个是曼宁糙率系数n。天然河道的n在 0.025 到 0.05 之间城市下垫面则在 0.015 到 0.12 之间。用 0.01 做初始化虽然跑得快但流速偏大、锋面推进过快最后实测对不上问题不在格式而在底摩擦没标定。第三个是时间步长限制中的底摩擦源项。显式格式里底摩擦项在浅水区会让局部特征波速变大若不对步长做限制在干湿边界附近极易振荡。所以在每个推进步里我都会加入底摩擦源项的特征频率限制if (h ds) then lambda_friction 2.0 * g * n2 * sqrt(u*u v*v) / (h**(4.0/3.0)) dt_limit min(dt_limit, 0.5 / max(lambda_friction, tiny(1.0))) end ifn2是曼宁系数的平方h是水深u、v是速度分量。这段逻辑的本质是把底摩擦看作一个隐式源项避免它在浅水区主导时间步长。放入主推进循环后矩形水槽算例的稳定步长往往可以提升 20% 以上。5. 把 typhon-solver 接入自己的网格与数据流5.1 网格格式转换技巧typhon-solver 的非结构网格输入往往要求“节点表 单元表 边表”三个数组。第三方网格生成工具Gmsh、Triangle、ANSA输出的格式五花八门我一般用一个小 Python 脚本做预处理把单元邻接关系重新编号为边表def build_edges(elems): edge_dict {} edge_list [] for cell, verts in enumerate(elems): for i in range(len(verts)): v1, v2 verts[i], verts[(i1) % len(verts)] key (min(v1, v2), max(v1, v2)) if key not in edge_dict: edge_dict[key] len(edge_list) edge_list.append([key, cell]) else: edge_list[edge_dict[key]].append(cell) return edge_list这段函数将每个内部边记录两个相邻单元编号边界边则只有一个单元编号。后续 Fortran 程序只要读取二维数组即可完成边表构建省去了在编译后的程序里做哈希查找的开销。需要注意的是坐标精度浮点坐标建议用双精度写出避免节点重合导致边表错乱。5.2 用 Gauge 点验证水位过程线输出场文件的验证成本太高我通常会在算例里埋 3 到 5 个 Gauge 点把每个时间步的水深和流速增量写入 CSV 文件。Gauge 点的时序数据可以直接和实测水位、解析解或商业软件结果对比判断收敛性和参数标定是否合理。输出频率不要太密建议至少每 5 个时间步写一次否则 I/O 会把整个计算拖慢三倍以上。5.3 一个快速自查技巧每次改动源码或重编后先跑一个静水池算例初始水深常数、无入流无出流、无地形坡度。程序应保持水位不变误差控制在机器精度量级。若静水池水位漂移超过 1e-6 m问题基本指向通量梯度与地形源项不平衡。把静水池测试通过后再跑真实地形省掉大半排错时间。新 mesh 上手时同时盯住启动后的前 100 步水位极值若持续振荡且振幅不衰减先把cfl减半若仍不收敛再去检查边界条件里的特征线方向是否配对。本文还有配套的精品资源点击获取