简介typhon-solver-0.3.0 是一套基于 Fortran 的二维浅水方程求解源码包采用非结构网格上的有限体积法面向从事洪水模拟、海洋动力学、河流与海岸水流计算的科研人员和工程师也适合数值方法初学者作为入门范本。压缩包共 378 个文件主体为 318 个 f90 源文件另有 32 个 make 构建脚本和若干辅助文件整体仅 766KB便于下载、阅读与改造。源码完整涵盖了非结构三角形网格生成、高阶数值积分、Runge-Kutta 时间推进、广义拉格朗日乘子法GCL以及多种边界条件处理并附有山洪暴发、河流流动等算例内容预览中的 interpol、mgrid、cgns 模块进一步展示了插值、网格生成及 CGNS 数据接口的具体实现。读者可通过研读这些代码理解二维浅水方程的离散求解全过程并学习 Fortran 科学计算程序设计技巧。目前已有 504 人学习下载资源小巧且结构清晰具有较高的参考价值。1. 从 typhon-solver 0.3.0 源码包看 Fortran 二维浅水方程求解器的定位拿到typhon-solver-0.3.0-sources.tar.gz多半是接手了洪水演进或溃坝分析计算域是不规则河道或城区网格只能切成三角形算力又撑不起三维模型此时 Fortran 写的二维浅水方程非结构网格求解器就是最现实的选择。这类源码包不带图形界面也不替你准备网格。能不能跑出结果取决于四件事方程怎么离散、Fortran 编译环境怎么搭、非结构网格文件怎么生成、时间步长与干湿参数怎么设。下面按这条路径推进适合没跑通正在排错的工程师也适合想看清这类模型内部组织方式的从业者。文中命令和代码取这类求解器的通用做法拿到你的包后对照源码调整即可。2. 二维浅水方程在非结构网格上的有限体积离散守恒形式与通量计算2.1 先把控制方程写成守恒形式二维浅水方程守恒形式的出发点是水深和动量。守恒向量取 U (h, hu, hv)^T其中 h 是水深u、v 是 x、y 向的垂线平均流速。x 向通量 F (hu, hu² gh²/2, huv)^Ty 向通量 G (hv, huv, hv² gh²/2)^T。源项 S (0, −gh·∂z_b/∂x − τ_bx/ρ, −gh·∂z_b/∂y − τ_by/ρ)^Tz_b 是底床高程τ_b 是床面切应力通常用曼宁公式闭合τ_bx/ρ g·n²·u·√(u²v²) / h^(1/3)n 是曼宁糙率系数。把方程写成 ∂U/∂t ∂F/∂x ∂G/∂y S 之后最关键的认识是这个系统在结构上和可压缩 Euler 方程完全同构gh²/2 对应压力项√(gh) 对应声速。这意味着 Godunov 型格式、间断捕捉、近似黎曼解这套工具可以原样搬过来typhon-solver 这类 Fortran 求解器里常见的 hllc、roe 关键字就是这么来的。阅读源码前还要先建立一条底线静水条件h z_b constu v 0下压力梯度和底坡源项必须严格抵消。这个平衡关系决定了底坡项怎么离散也是后面排错时最先怀疑的对象。2.2 非结构网格的拓扑组织节点、单元、边三张表非结构网格没有 (i,j) 索引可依赖拓扑在代码里通常用三张表表达节点坐标表 node_xy、单元-节点表 cell_node、边-左右单元表 edge_lr。二维浅水几乎都用三角形网格因为 Delaunay 剖分成熟、能贴复杂河岸局部加密方便四边形剖分在跨河断面附近容易出畸变单元除非代码明确支持不建议首选。最小数据结构写成 Fortran 是这样module mesh_mod ! 三角形非结构网格的最小拓扑 integer, parameter :: dp kind(1.0d0) integer :: ncell, nedge, nnode real(dp), allocatable :: node_xy(:,:) ! (2,nnode) 节点坐标 integer, allocatable :: cell_node(:,:) ! (3,ncell) 单元三顶点 integer, allocatable :: edge_lr(:,:) ! (2,nedge) 左/右单元编号 real(dp), allocatable :: cell_area(:) ! 单元面积 real(dp), allocatable :: cell_xy(:,:) ! (2,ncell) 单元中心 contains subroutine mesh_read(node_file, elem_file) character(*), intent(in) :: node_file, elem_file integer :: i open(10, filenode_file, statusold) read(10, *) nnode ! 首行写节点总数 allocate(node_xy(2,nnode)) do i 1, nnode read(10, *) node_xy(1,i), node_xy(2,i) end do close(10) ! elem 文件同理之后算面积和中心 end subroutine mesh_read end module mesh_modedge_lr 是有限体积格式的中枢任一条边左单元 L、右单元 R数值通量定义为从 L 穿过边流入 R。边界边把其中一侧置 0 或负值并在另一个数组里记录边界类型固壁、开边界。cell_area 在守恒量更新时做分母如果某单元面积为零或负值头几步就会输出 NaN——这是后面专门要排查的网格质量问题。单元中心 cell_xy 在 MUSCL 重构时用来算距离和梯度不能省略。2.3 边法向通量HLL/HLLC 选型与限制器对每条边取法向单位向量 n (nx, ny)单元中心量通过重构外推到边中点得到左右状态 q_L、q_R再交给近似黎曼解计算法向通量。波速估计最常用 Davis 形式S_L min(u_L − c_L, u_R − c_R)S_R max(u_L c_L, u_R c_R)其中 u_L u_L·nx v_L·ny 是法向速度c_L √(g·h_L) 是浅水波速。得到三波速后中间波速由守恒条件确定。HLLC 的中间波速 S_* 在代码里常见这样写! HLLC 中间波速来自跨中间波通量的一致性条件 Sstar (SR*hR*unR - SL*hL*unL - 0.5d0*g*(hR*hR - hL*hL)) / (SR*hR - SL*hL 1.0d-30)说明分子里 0.5g(hR² − hL²) 是两个单元压力项之差分母上的 1.0d-30 防止干单元下 SL·hL 与 SR·hR 同时为零造成除零。注意这里必须用双精度字面量如果包里有单精度编译宏1.0d-30 会被截断成 0除零后直接 NaN。三种格式怎么选参考下表。保守经验调试期用 HLL 最稳正式方案比选用 HLLC。格式波数中间态浅水适用性HLL2单一常数态最鲁棒接触间断被抹平干湿界面不易出负水深HLLC3两个中间态精度更高浅水主流选择需配合干湿阈值Roe线性化无中间态间断分辨率好需要熵修正干底时容易负水深提示拿到源码包后不用通读全部代码先 grep 通量子程序里的 hllc、roe、minmod 字符串几分钟就能定位格式、限制器和干湿处理的位置。3. 编译 typhon-solver 0.3.0Fortran 环境安装与 sources.tar.gz 构建流程3.1 Fortran 编译器安装与版本验证Linux 下最常用的 Fortran 编译器是 gfortran安装方式# Debian / Ubuntu sudo apt update sudo apt install -y gfortran # CentOS / RHEL sudo yum install -y gcc-gfortran # 验证版本 gfortran --version装完别急着编译先验证默认精度和标准符合性。老代码常见 real(kind8)、REAL*8、implicit none 混用新版本 gfortran 对标准符合性检查更严00 年代代码经常报一堆警告。用一行命令确认工具链可用echo print *, ok, selected_real_kind(15) | gfortran -x f95 - -o /tmp/fcheck /tmp/fcheck如果包带了 MPI 并行需要额外装 mpich 或 openmpi并通过 mpifort 包装器编译sudo apt install -y openmpi-bin libopenmpi-dev export FCmpifort mpifort --version这里的选型原则串行调试用 gfortran跑正式算例再切 MPI 版本不要在调试期就让并行环境干扰排查。安装完编译器后建议把--version输出里的小版本号记下来后面遇到浮点行为异常时编译器版本和优化级别是最先要对比的变量。3.2 解包、看构建脚本并完成编译tar xzf typhon-solver-0.3.0-sources.tar.gz cd typhon-solver-0.3.0 find . -maxdepth 2 -type f | sort | head -40解包后先看 README 和构建脚本比看源码省时间。这类 Fortran 项目最常见的构建方式有两种纯 Makefile 或 CMake。Makefile 项目的标准做法make FCgfortran FFLAGS-O2 -g -fcheckbounds -Wall -j4CMake 项目的标准做法cmake -B build -DCMAKE_Fortran_COMPILERgfortran \ -DCMAKE_BUILD_TYPERelease -DCMAKE_Fortran_FLAGS-O2 cmake --build build -j4FFLAGS 里每个选项都有明确的用途和生命周期参数作用何时开/关-O2常规优化级别正式计算至少 -O2-g生成调试信息需要 gdb 或 backtrace 时开-fcheckbounds数组越界检查调试期开跑大网格关掉-Wall全部警告全程开警告多说明隐式类型或未初始化-ffree-line-length-none允许任意行长老代码常见长行默认 132 列会报错编译报错最常看见三类。第一类是undefined reference to ... gfortran_*多半是链接器用了 gcc 而不是 gfortranMakefile 里 FC 和 LINKER 不一致。第二类是Type mismatch in argument典型原因是某个子例程没有 use 对应模块触发隐式类型规则顺着调用链到子例程声明处就能看见。第三类是Syntax error in OPEN statement常见于 STATUSNEW 在文件已存在时的处理新编译器直接报错。3.3 编译产物自检与输入方式确认编译完成后产物可能是单个可执行文件也可能是 bin/ 下多个工具网格检查器和求解器分开。先用无参数方式跑一次看帮助./typhon-solver --help 21 | head -20如果 README 太长直接提取 quick start 段落grep -A 20 -i quick start\|usage README这一步的核心目的是确认可执行文件的输入方式是命令行传配置还是标准输入读 namelist。多数 Fortran 求解器两者都支持初跑建议用 namelist 方式查参数比扒命令行选项直观得多下一章展开配置结构。4. 准备非结构网格与输入配置跑通第一个二维浅水算例4.1 用 Gmsh 生成三角形网格的完整流程第一个算例建议从矩形域加一个突起的溃坝开始网格和边界都好控制。Gmsh 的 .geo 文件// bump.geo40 m x 10 m 矩形河道最大单元尺度 0.5 m SetFactory(OpenCASCADE); Rectangle(1) {0, 0, 0, 40, 10}; Mesh.CharacteristicLengthMax 0.5; Mesh.Algorithm 6; // Delaunay 三角化 Mesh 2;生成命令gmsh bump.geo -2 -o bump.msh。说明Mesh.Algorithm 6 对应 Delaunay二维浅水几乎都用它生成的三角形质量均匀复杂河道地形后期改用背景场控制单元尺度这里先用均匀网格验证流程。很多 Fortran 求解器不直接读 msh而是要 node/elem 两个纯文本文件。转换脚本用 meshio 最省事import meshio import numpy as np mesh meshio.read(bump.msh) np.savetxt(bump.node, mesh.points[:, :2], fmt%.10f) tri [c.data for c in mesh.cells if c.type triangle][0] np.savetxt(bump.elem, tri 1, fmt%d) # 节点编号从 0 改 1提示下标 1 是非结构网格格式转换里最常见的坑。Gmsh 的 msh 节点从 0 计数Fortran 数组默认从 1 计数漏掉这一步第一个单元会引用 0 号节点读进去就是越界或负面积。转换前还要确认一件事文件首行要不要写节点数、单元数。看源码里 mesh_read 的 read 语句前面读几个量就知道这决定了生成文件的头部格式。4.2 namelist 配置结构与参数选择Fortran 求解器最常见的输入方式是 namelist跟语言内建 I/O 天然配合。配置通常长这样mesh node_file bump.node elem_file bump.elem bc_file bump.bc / physics grav 9.81 manning 0.025 t_end 30.0 / scheme riemann hllc limiter minmod cfl 0.8 theta_wet 1.0e-3 / output interval 100 format vtk output_dir out /参数说明表参数典型范围备注cfl0.5–0.9显式格式硬约束超过 1.0 几乎必然发散theta_wet1e-3–1e-5水深低于此值视为干单元manning0.01–0.06混凝土 0.013天然河道 0.03–0.05grav9.81不要顺手写 9.8验证算例对不上解析解注意参数名每个包都不一样。拿到包后先在源码目录执行grep -rn namelist .找到声明块里的变量名再对齐 nml 文件否则 namelist 读到未知项会直接跳过参数没生效但你完全不知道。4.3 运行算例并确认推进节奏mkdir -p out ./typhon-solver config.nml run.log 21 tail -f run.log运行日志重点看两类信息时间推进行和守恒量报告。正常情况每步打印步号、时间、水深范围类似step 100 t 2.53123 hmin0.0000 hmax2.8471 mass1.234e05hmax 应该从初始值缓慢变化而不是跳变mass 在没有进出流时基本不动。输出是 VTK 就进 ParaView 看水深云图二维浅水算例首先要确认的只有一件事水有没有从该走的方向流过去。5. 稳定性控制与排错二维浅水算例的 CFL、干湿界面与 NaN 定位5.1 时间步长的 CFL 约束与全局归约显式有限体积格式的时间步长受网格尺度约束信息最多在一个步长内跨越一个单元。非结构网格下每个单元的局部步长是 dt_i CFL·L_i / (|u_i| c_i)L_i 的取法影响很大保守用单元最小边长想提高效率用内切圆直径。当网格单元大小悬殊时全局步长被最小单元拖累这就是为什么网格质量直接决定算例能不能在合理时间内跑完。全局步长归约在 Fortran 里常见这样写! 全局时间步长归约取所有单元局部步长的最小值 dt_global huge(1.0d0) !$omp parallel do reduction(min:dt_global) do i 1, ncell c sqrt(g * max(h(i), 0.0d0)) ! 防负水深 dt_global min(dt_global, cfl * edge_min(i) / (abs(u(i)) c 1.0d-12)) end do !$omp end parallel do说明max(h, 0.0d0)防止对负水深开根号1.0d-12是干单元下的防除零保护量级要比正常波速小几个数量级才不会干扰步长用reduction(min:...)而不是 critical 区块OpenMP 下归约效率高得多这是很容易被忽略的并行性能点。5.2 底坡源项平衡与干湿阈值联动底坡源项是这类求解器里最容易被写错的代码段。静水条件下通量里的压力梯度和底坡项必须严格抵消。验证手段是放一个倾斜底床、初始完全静水的算例跑上千步如果最大流速超过 1e-10 量级说明源项平衡有问题常见修法是 hydrostatic reconstruction按边两端底高程重构水位而不是直接插水深。干湿阈值 theta_wet 要和网格尺度联动。设太大真实的浅水前沿会被削平水位提前出界设太小干单元先出现轻微负水深随后一步被波速计算放大成 NaN。经验取值是当地水深量级的 1e-4 到 1e-3并配合h max(h, 0.0d0)的抹零逻辑一起用。调试期内宁可让阈值偏大先拿到稳定解再往下收紧。5.3 按发生步数定位 NaN 来源NaN 的出现位置本身就是诊断信息按发生时间分类现象最常见原因先查哪里第 1 步就 NaN单元面积零或负、初始场负水深网格面积检查、初始水位赋值固定步数后才 NaNCFL 偏大或波速估计在干湿界面失稳把 cfl 降到 0.4 复跑对比个别单元局部 NaN网格畸变细长三角形、面积悬殊检查该单元最小角和面积比dt 骤减到接近 0局部流速异常大水深被抹成负在 h 0 处打断点调试通用做法编译保留-g -fcheckbounds在守恒量更新之后加一段哨兵检查do i 1, ncell if (h(i) / h(i)) then write(*,*) NaN at cell, i, time , t stop end if end do提示h(i) / h(i)是 Fortran 里判断 NaN 的标准写法IEEE 754 规定 NaN 不等于任何数包括它自己。找到出问题的单元编号后把它的三个节点坐标、相邻单元面积打出来十有八九指向网格质量问题而不是格式问题。这个定位习惯能省掉大量瞎调参数的时间。6. 结果验证与性能优化把算例从能跑推到可信6.1 质量守恒与解析解双校验质量守恒是最低门槛的验证。从 VTK 序列读水深和面积算总水量import meshio, numpy as np, glob files sorted(glob.glob(out/*.vtu)) for f in files: m meshio.read(f) h m.cell_data[h][0] pts, tri m.points, m.cells[0].data a np.abs(np.cross(pts[tri[:,1]] - pts[tri[:,0]], pts[tri[:,2]] - pts[tri[:,0]])) / 2.0 print(f.split(/)[-1], float(h a))守恒格式下总水量偏差应该在机器精度量级实际偏差主要来自边界进出流。如果偏差随时间线性增大优先怀疑边界通量实现而非格式本身。接着做解析解对照一维溃坝用 Ritter 解二维圆形溃坝验证对称性。做法是把水深沿半径提出来和精确解画在一起网格加密后最大 L2 误差应接近二阶这是 MUSCL-HLLC 组合的理论预期。6.2 编译选项与 OpenMP 并行实操热点是边循环的通量计算。常见做法是先按边算通量存数组再按单元做守恒量更新这样 OpenMP 并行边循环没有写冲突export OMP_NUM_THREADS8 ./typhon-solver config.nml编译用-O3 -marchnative把本机指令集用满但始终不要开-ffast-math——它会破坏 NaN 判断语义并降低精度浅水格式的干湿判定和哨兵检查全部依赖 IEEE 语义。开-fopenmp后注意时间步归约必须写 reduction 而不是 critical否则 8 线程可能比单线程还慢。最后一个实用验证技巧同一套网格上把 cfl 从 0.8 降到 0.4 各跑一版两张水深等值线如果完全重合说明时间离散误差不是主导项可以放心用大步长做批量方案比选如果两张图明显分开先查限制器和重构精度再考虑压缩步长。本文还有配套的精品资源点击获取