ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

Basilisk气泡模拟实战:C语言宏与Shell容器化

Basilisk气泡模拟实战:C语言宏与Shell容器化 简介一份围绕 Basilisk 开源数值模拟框架整理的博士课题学习资料面向流体力学、地球物理等方向、需要借助 C 语言与 Shell 脚本完成仿真研究的科研人员也适合从偏微分方程数值求解入门到进阶的个人学习者。资料从有限体积法和谱方法等基础入手串联科学计算、编译器、调试、版本控制、并行计算与可视化等环节可帮助梳理从方程离散化到自动化模拟运行的完整链路。压缩包共 14 个文件约 508KB以 C 源文件、头文件、Shell 脚本和 Makefile 为主另含 README 说明文档、LICENSE 与 .gitignore 配置。C 文件与头文件用于定义求解器并输出 VTU 格式结果Shell 脚本负责容器环境启动与参数扫描等任务自动化Makefile 简化构建流程。目前已有 211 人浏览学习。通过示例代码和配置脚本可以借鉴环境搭建、编译调试、并行计算和结果分析的一套完整思路对开展气泡轴对称/二维模拟等博士研究具有直接参考价值内容来自网络分享使用前请自行核对版权与可用性。1. Basilisk模拟与我的博士研究C_Shell方案怎么用博士课题做到多相流阶段我拿到一个名为“与我的博士相关的Basilisk模拟_C_Shell_下载”的压缩包里面不是现成的可视化软件而是一套Basilisk源码工程两个气泡求解器文件、一个Makefile、两个容器启动脚本以及几个后处理头文件。Basilisk是开源CFD框架用C语言编写配合Shell脚本完成编译、运行和结果整理。它的真正门槛在于C语言宏扩展和事件机制第一次打开代码会感觉不像标准C但只要理解了分层逻辑就能很快改出自己需要的模拟场景。这篇内容面向博士研究生、科研人员以及需要用Basilisk做两相流的工程师从框架原理讲到容器化运行和调参验证中间会给出可以直接照抄的命令和参数表。2. Basilisk框架与C语言实现树状网格、VOF和事件循环2.1 C语言宏Basilisk的“伪面向对象”机制Basilisk的核心不是传统CFD程序那种类封装而是通过C预处理宏在编译期生成数据结构。以压缩包中的bubble_2D.c为例代码里会出现scalar f[];这样的声明这里的f不是定长数组而是一个标量场。宏在编译展开后变成带网格单元信息的结构体通过foreach循环遍历网格单元。这种设计让求解器代码接近纯C性能同时允许用户在问题文件里像写脚本一样定义物理场。这里有一段常见的Basilisk风格C代码示意了如何在一个简单的扩散步中操作场变量#include grid/octree.h #include navier-stokes/centered.h scalar T[]; face vector flux[]; event diffusion_step (i) { foreach() T[] dt * (flux.x[] flux.y[] - flux.x[1,0] - flux.y[0,1]) / sq(Delta); boundary ({T}); }代码中flux.x[]表示单元左侧面的通量flux.x[1,0]表示右侧面的通量flux.y[]和flux.y[0,1]则对应下上两个面。Delta是当前网格单元的尺寸dt是时间步长sq()是Basilisk提供的平方宏。foreach()遍历当前层级的所有网格单元boundary({T})负责在自适应网格交界处同步字段值。Basilisk的索引语法看起来像标准C数组实际上通过宏映射到了树状网格单元理解这一点后再去读bubble_axis_symmetric.c中的VOF输运代码就不会发怵。2.2 事件循环模拟节奏的控制核心Basilisk程序通常由main()设置物理参数和初始网格然后调用run()启动事件调度。事件用event关键字定义括号里写触发条件。为什么不用普通的while循环因为科研模拟往往需要在指定时间点输出、自适应加密、改变物理规则拆成事件可以让这些逻辑互不干扰。典型的初始化与日志事件如下event init (t 0) { fraction (f, sq(x) sq(y) - sq(0.1)); } event logfile (i) { fprintf (stderr, %g %g %d\n, t, dt, grid-tn); }event init在模拟开始前执行fraction函数根据隐式方程sq(x)sq(y)-sq(0.1)0标记初始气泡区域。event logfile在每一个迭代步后执行i表示“每次迭代都触发”。grid-tn是当前网格总单元数能实时反映自适应加密是否工作。事件机制把时间推进、网格自适应和文件输出解耦这也是同一个求解器内核可以套用在完全不同的物理问题上的原因。2.3 VOF方法与自适应网格的组合气泡模拟属于典型的两相流Basilisk的two-phase.h模块使用几何VOF方法传输体积分数f。界面处的密度和黏度由f平均表面张力通过连续表面力模型加入动量方程。自适应加密在界面附近提高分辨率可以大幅降低计算量。加密逻辑通常放在一个独立事件里event adapt (i) { adapt_wavelet ((scalar *){f, u.x, u.y}, (double[]){1e-2, 1e-2, 1e-2}, 10, 0); }adapt_wavelet是Basilisk的小波误差控制函数第一个参数指定要加密的字段第二个参数指定每个字段的绝对容差第三个参数是最大网格层级MAX_LEVEL第四个参数是最小层级。容差越小网格越细最大层级每增加1理论网格量增加4倍二维。这里有一个常见误用把所有字段容差设成相同值。速度场和体积分数的尺度差异很大实际调试时建议分开设置例如f用1e-3速度场用1e-2否则会出现界面附近已经很细但速度梯度大区域仍然捕获不足的情况。树状网格上的通量和梯度插值比结构化网格复杂Basilisk用face vector、gradient等宏封装底层细节。对使用者来说主要收益是可以用很少的代码完成“界面加密、远处粗化”的网格策略这也是Basilisk在微流体和气泡模拟中流行的重要原因。3. Shell脚本与容器化createContainer.sh和startContainer.sh的使用逻辑3.1 为什么科研项目要把Basilisk装进容器Basilisk在编译时依赖特定版本的GCC、OpenMPI、GSL以及ffmpeg。不同Linux发行版自带的工具链版本不同很容易出现编译通过但运行崩溃的情况。我在不同服务器上遇到过几次Ubuntu 22.04上编译好的可执行文件放到CentOS 7上直接报缺少libgomp符号。压缩包里的createContainer.sh和startContainer.sh就是为这个问题设计的它们把工具链和运行环境固定在容器里宿主机的系统库版本不再影响模拟结果。常见做法是先在宿主机上准备Docker或Podman环境然后由createContainer.sh构建一个包含固定编译工具的镜像startContainer.sh负责创建并进入运行容器同时把当前工作目录挂载进去。3.2 createContainer.sh构建可复现的编译环境从脚本命名和博士项目习惯看createContainer.sh负责生成镜像。一个典型实现是这样的#!/bin/bash set -e IMAGEbasilisk/phd:1.0 docker build -t $IMAGE -f Dockerfile.basilisk . echo image $IMAGE readyset -e是Shell脚本里的保险丝任何一条命令返回值非零都会立即终止脚本避免用一个损坏的镜像继续往下跑。Dockerfile.basilisk里一般会基于Debian或Ubuntu镜像安装build-essential、git、gnuplot、libgsl-dev等依赖然后克隆Basilisk源码并设置环境变量BASILISK。这样构建出来的镜像就是可重复的编译环境。这里有个细节值得注意如果docker build执行到一半因为网络原因失败set -e会阻止脚本继续但已经下载的镜像层会留在本地。再次执行构建时Docker会使用缓存速度会快很多不是bug。3.3 startContainer.sh挂载、权限和并行参数进入模拟阶段时startContainer.sh启动交互式容器。从使用场景看它需要把源码目录挂载进去同时避免在容器里产生root权限的文件。一个可用模板是#!/bin/bash set -e IMAGEbasilisk/phd:1.0 WORKDIR/sim docker run --rm -it \ --name basilisk_run \ -v $(pwd):$WORKDIR \ -w $WORKDIR \ -u $(id -u):$(id -g) \ -e OMP_NUM_THREADS${OMP_NUM_THREADS:-4} \ $IMAGE bash几个关键点--rm保证退出后容器被删除不占用磁盘-v $(pwd):$WORKDIR把当前目录挂载为容器工作目录模拟产生的结果文件都落在宿主机方便后续用Paraview查看-u把当前UID和GID写入容器避免生成root权限的结果文件-e OMP_NUM_THREADS控制容器内使用的OpenMP线程数。进入容器后就可以按通常流程执行make ./bubble_2D。如果容器内没有编译过源码第一次执行make会比较慢因为Basilisk会把所有依赖头文件都扫描一遍。后续再编译时只有改动的C文件会被重编。3.4 Shell脚本参数速查下面这张表总结了脚本里最常调整的参数我平时会贴在项目目录下Shell语句作用注意点set -e命令失败即退出管道场景下要用set -o pipefail配合docker build -f Dockerfile.basilisk按指定文件构建镜像Dockerfile文件名若非默认路径必须指定docker run --rm容器退出后自动删除配合-d后台运行时不适用-v $(pwd):/sim挂载当前目录路径含空格时务必加引号-u $(id -u):$(id -g)避免生成root文件Windows下要改用Docker Desktop的权限映射-e OMP_NUM_THREADS设置OpenMP线程数线程数应小于等于物理核心数这张表实际使用时会反复查。一个常见问题是在容器中运行mpirun还需要把宿主机的高速网络接口映射进去通常需要--network host。但博士课题的模拟单机OpenMP已经够用优先不要上MPI因为Basilisk的MPI编译要比OpenMP更容易踩坑而且泡模拟的单节点算力通常足够支撑十万到百万级网格。4. 气泡模拟实战bubble_axis_symmetric.c与bubble_2D.c编译调参4.1 轴对称与二维模型的选型差异压缩包里的bubble_axis_symmetric.c和bubble_2D.c代表两种不同维度的模拟方式。bubble_2D.c在二维矩形域里模拟圆形气泡的上升计算量小适合做参数扫描和数值实验bubble_axis_symmetric.c通过axi.h把二维域解释为轴对称坐标得到的流场更接近三维球形气泡。选型时的一个基本原则如果只关心单个气泡的终端速度和形状优先用轴对称模型它兼顾了二维的计算效率和三维物理特征如果研究气泡对之间的相互作用轴对称就无法描述非轴对称扰动必须用二维或三维全模型。bubble_axis_symmetric.c里的重力方向通常沿轴向半径方向由y坐标表示所以初始气泡形状的定义要写成sq(x)sq(y)而不是sq(x)sq(z)。4.2 Makefile与编译命令Basilisk使用自带Makefile体系。压缩包中的Makefile一般会写成BASILISK ? $(HOME)/basilisk include $(BASILISK)/Makefile.defs EXECS bubble_2D bubble_axis_symmetric bubble_2D: bubble_2D.c bubble_axis_symmetric: bubble_axis_symmetric.c include $(BASILISK)/Makefile.rules这里的BASILISK环境变量必须指向Basilisk源码目录否则使用默认的$(HOME)/basilisk。Makefile.defs负责探测编译器和库Makefile.rules补充通用编译规则。编译命令很直接make bubble_axis_symmetric OMP_NUM_THREADS4 ./bubble_axis_symmetric log 21第一条命令生成可执行文件。第二条命令用4线程运行并把标准和错误输出都重定向到log。Basilisk的进度信息默认写到stderr所以如果你只重定向标准输出tail -f log里会什么都看不到。编译时如果出现与fractions_output.h相关的错误通常是因为当前目录下缺少Basilisk头文件检查BASILISK路径是否设置正确。提示Basilisk的标准输出默认不记录进度必须把stderr一起重定向否则tail -f log里永远是空的。4.3 参数表与核心调参逻辑下列参数在两个C文件中基本都能找到只是初始值不同参数作用调节方向size()计算域边长至少为气泡直径的6倍否则边界效应明显init_grid()初始网格数通常取1 6配合后续自适应加密rho1, rho2液体/气相密度密度比过大会导致时间步长急剧变小mu1, mu2液体/气相动力黏度注意Basilisk默认是无量纲方程f.sigma表面张力系数减小会抑制界面变形为零则无界面张力MAX_LEVEL最大网格层级每加1最小网格尺寸减半CFL时间步长安全系数默认约0.5发散时降到0.2G重力加速度控制无量纲佛洛德数调参有一个常见误区直接改MAX_LEVEL就认为会“更准确”。网格加密后会解析出更小的涡和界面结构但如果初始扰动或边界条件没有相应变化结果可能更发散。我的做法是先固定其他参数只把MAX_LEVEL从6逐步升到10比较气泡质心速度和界面最大曲率的收敛情况。CFL参数也需要配套调整网格变细后同样的物理时间步需要更多迭代步如果CFL取得过大显式格式会在计算下一时间步时产生nan。4.4 运行时的坑与定位手段实际跑气泡模拟最容易遇到三种情况。第一种是程序启动后很快输出nan大概率是初始f.sigma偏小或CFL偏大表面张力模型无法稳定支撑界面。可以先检查log中最后一个完整的dt值再把CFL降到0.2重跑。第二种是网格一直加密到最大层级后程序变慢但不报错这说明adapt_wavelet的容差设得过小可以放宽速度场容差例如把u.x, u.y的容差从1e-3放宽到1e-2。第三种是make编译时提示fatal error: common.h: No such file or directory这通常是因为BASILISK环境变量没有在当前Shell会话里生效执行export BASILISK/path/to/basilisk make clean make bubble_axis_symmetricmake clean必须执行否则已经生成的依赖文件会继续引用旧的包含路径即使环境变量修正了编译仍会失败。5. 进阶VTU输出、并行计算与网格无关性验证5.1 把模拟结果导出为Paraview可读的VTU文件压缩包中的output_vtu_foreach.h就是用于输出VTU的。Basilisk官方库里既有output_vtu.h也有output_vtu_foreach.h后者更适合在OpenMP并行下写文件因为它把输出动作放在遍历循环内避免多个线程同时写同一个文件描述符。在C文件里加上一个事件即可#include output_vtu_foreach.h event vtu_output (t 0.1) { char name[80]; static int fnum 0; sprintf (name, snap-%03d.vtu, fnum); FILE *fp fopen (name, w); output_vtu_foreach ((scalar *){f, u.x, u.y}, (char *){f, ux, uy}, fp); fclose (fp); }t 0.1表示每增加0.1时间单位输出一个VTU文件fnum保证文件名不重复。输出的文件可以直接用Paraview打开通过Threshold过滤器提取气液界面。注意这里输出的是Basilisk内部的坐标一般是无量纲的与实验对照时必须在后处理阶段乘以特征长度。5.2 并行计算的线程数选择在startContainer.sh里通过环境变量设置OMP_NUM_THREADS但Basilisk运行时读取的就是该变量。线程数不是越多越好自适应网格的负载分布并不均匀线程过多会导致内存带宽竞争。单机的经验阈值是8线程以内加速比接近线性超过后收益明显下降。如果要在多节点上跑就必须启用MPI重新编译Basilisk这不是单机项目优先考虑的事。5.3 网格无关性验证的快速方法最后给一个可用的验证思路用脚本循环改写源码中的最大网格级别跑多个分辨率然后比较气泡质心高度随时间的变化。可以这样组织for level in 8 9 10; do sed -i s/#define MAX_LEVEL [0-9]*/#define MAX_LEVEL $level/ bubble_2D.c make bubble_2D /dev/null 21 OMP_NUM_THREADS4 ./bubble_2D 2level_$level.log donesed替换源码中MAX_LEVEL的值前提是文件里已经有#define MAX_LEVEL这一行。如果没有就在main()函数前加一行。跑完后用awk提取三个日志中气泡质心速度的稳态值检查相邻两个级别的相对误差。如果相邻两个分辨率的稳态速度差异小于1%就直接用较细的那个级别出结果省下的算力足够再多跑一组参数扫描。本文还有配套的精品资源点击获取
返回列表