OpenFOAM中的二维与三维自适应网格负载均衡
摘要
本文针对OpenFOAM（一个成熟的C++计算流体动力学库）提出了一种面向对象的二维和三维负载均衡自适应网格细化方法。这种高性能计算技术对于采用域分解的分布式内存并行计算机架构而言，在处理移动关注区域时是实现计算效率的必要手段。移动关注区域会在模拟过程中动态变形并跨越计算域迁移，这要求对解的特征保持高空间分辨率。我们详细阐述了软件设计与代码结构，重点介绍了在工程应用软件可用性方面取得的成果，以及为实现高稳定性和运行时性能所必需的多项技术改进。
1 引言
科学与工程领域中瞬态输运过程的并行计算模拟常常面临一个重大挑战：如何在计算过程中高效利用计算资源，同时始终保持数值解的高精度。输运过程往往涉及动态变化的关注区域，这些区域会随时间推移在计算域内迁移并发生形变。为降低数值误差，需要在这些动态区域内完全解析细微的瞬态解特征。工程应用中的典型实例包括火焰锋面、爆震激波或流体界面。
当前，计算流体力学（CFD）软件通常依赖分布式内存并行计算机架构来处理所有中大规模计算任务。现代CFD仿真软件要求其基础核心—有限体积法（FVM）必须支持具有任意拓扑结构的动态非结构化计算网格，以应对形状多变的复杂求解域。在此类需求下，业界普遍采用计算网格上的同位（伪交错）变量排列方案：输运变量值存储于单元中心（体积场），而单元面中心则保存通量信息及插值单元值（表面场）。并行化处理几乎完全通过区域分解模式实现—即将网格单元面的大型循环任务分割后交由多个处理器核心（下文简称处理器）协同处理。
目前已有多个商业和学术CFD求解器库具备自适应网格细化(AMR)和动态负载均衡(DLB)功能。在已提及的DLB实现方案中，本研究具有重要意义，因为这OpenFOAM首个稳定实现版本。利用本文提供的AMR与DLB框架开展的仿真研究已在文献和中得到应用。
AMR与DLB的实现取决于底层数据结构和并行化策略。OpenFOAM（开放场运算与操作）是一个功能全面的开源C++库，专注于计算连续介质力学领域（包括CFD）。该库采用严谨高效的面向对象方法实现域分解并行计算，其模块化代码结构使其区别于其他开源CFD软件—线性求解器、插值方案、物理模型等模块均遵循策略设计模式实现。通过泛型编程实现可维护性：大量使用模板可在不同数据类型上执行相同操作而无需代码重复，结合运算符重载技术，能构建高度贴合连续介质建模数学语言（即偏微分方程）的顶层CFD代码。运行时类型选择(RTS)机制进一步增强了灵活性：遵循工厂设计模式，用户仅需修改配置文件（字典）中的条目，即可在软件运行时实例化类层次结构中的对象。这种接口与库使用之间的严格区分贯穿整个库设计，旨在隐藏实现细节以确保顶层代码的可用性和可读性。
在并行计算方面，通过作为并行通信封装器的软件接口层，将通信细节与顶层库使用隔离。该层提供标准接口，使得所有顶层代码的编写都无需考虑特定并行化要求，从而实现串行与并行执行使用相同代码。处理器间通信被建立为边界条件，采用零光晕层方法将每个网格单元唯一分配到不同处理器，避免在处理器边界处重复存储单元数据。
OpenFOAM的许多发展都是由社区驱动的，这极大地促进了该框架的成功。本研究正是通过修订并构建社区在二维和轴对称案例以及动态负载平衡(DLB)[16]方面的贡献，展现了这种集体努力的典范。本文将自适应网格加密(AMR)与DLB的改进整合为一个统一框架，并向社区开放。其优势在于能有效降低中小型工程问题的计算量—这类问题在现代台式机上即可求解，而AMR与DLB技术在此类场景中同样不可或缺。由于本文提出的AMR-DLB框架受益于众多OpenFOAM库，其并行性能无法脱离OpenFOAM环境独立测试。因此，详尽的并行效率与扩展性测试（尤其是面向更大规模众核架构的测试）不在本研究范围内。不过值得注意的是，Phuc等学者近期在静态加密网格上使用数十亿单元的实验已展现出令人鼓舞的结果。
为有效应对上述挑战，我们采用并显著增强了OpenFOAM C++库。本研究的贡献基于以下目标：
• 在OpenFOAM软件的面向对象架构中阐述自适应网格加密（AMR）与动态负载平衡（DLB）的核心要素
• 向社区提供可直接用于新版OpenFOAM（4.x及以上版本）多种现有求解器的实施方案
2 自适应网格细化
为实现仅在必要区域动态获取高精度，AMR（自适应网格加密）技术采用了h-自适应方法。目前OpenFOAM仅支持三维六面体网格的AMR功能。近期研究提出了适用于任意多面体网格的AMR扩展方案，该方案目前仅实现在foam-extend-4.1分支中尚未并入OpenFOAM主代码库，而文献的代码当前未公开获取。多面体网格细化会将单元分解为四面体，这会在二维和2.5维案例的恒定模式下产生寄生电流。因此，真正的二维及2.5维网格细化方案在许多应用场景中仍具有重要价值。
2.1 八叉树单元细化
图1 (a) 六面体单元的八叉树细化会生成12个内部面。表面场值（圆形符号）通过相邻主控面（方形符号）的数值进行插值计算，图中示例性展示了一个内部面的情况。(b) 六面体单元与棱柱单元的四叉树细化。该不变方向垂直于页面。
在三维情况下，六面体网格单元采用八叉树结构进行分割，每个单元关联一个所谓的细化层级。通过将母单元切割成八个子单元来实现细化，这些子单元共具有36个面，其中12个面位于母单元内部（见图1(a)）。与细化单元相邻的单元会从六面体转变为具有超过六个面的多面体单元。
OpenFOAM中的伪交错方法依赖于两个位置来定义网格相关场：体积场表示单元中心的值，而表面场表示面中心的值。当前OpenFOAM版本中，结合自适应网格细化（AMR）处理表面场时存在三个主要问题：细化与未细化单元值之间缺乏一致的插值方法、细化过程中新创建面的错误寻址，以及非通量表面场的符号翻转问题。下文将概述这些问题的解决方案。
映射是指父单元与细化子单元之间场的插值过程。单元中心体积场的映射以直接且保守的方式实现。在单元细化步骤中，子单元继承父单元的单元中心值；在单元反细化步骤中，父单元中心设置为子单元中心值的体积平均值。然而，细化步骤中表面场的映射更具挑战性和关键性，因为OpenFOAM中大多数求解器基于通量，而通量以表面场形式表示。
图2 UML图展示了新实现的动态负载均衡类、二维及2.5维网格细化类（以黄色标注），以及需要修改的既有类（以红色标注）。
细化面上的值设置为对应母面（称为主面）的值，主面是指被分割成四个新面的父单元面。但位于父单元内部的新面与任何主面都没有关联。本工作通过对四个相邻主面值进行算术平均，实现了对这些内部面更好的值插值，如图1(a)中一个内部面所示。表面场到新内部面的映射已在dynamicRefineFvMesh类中实现（见图2）。
修正了网格切割类hexRef8中可追溯至OpenFOAM 1.x版本的面初始化错误，该错误导致新创建内部面的寻址失效，使得整个模拟域中新面的值呈现看似随机的情况。
通量的符号取决于面法线方向，该方向由最多两个相邻单元共享。由于拓扑变化（在细化或粗化过程中），该方向可能会根据形式所有权发生改变，从而导致通量符号翻转。数值修正变得必要。符号翻转操作适用于所有类型的表面场，即使这些场可能并不代表通量。OpenFOAM库中多处存在的错误符号翻转问题已得到修正。
OpenFOAM中的不可压缩流求解器依赖于保守通量场，这对应于无散度速度场。由于表面场的非保守映射特性，每个自适应网格加密（AMR）步骤后都需要对通量场进行修正（见附录A）。通过采用改进的面寻址修正、增强的表面场映射以及符号翻转修正（参见提交记录1cef6f9），显著减少了确保无散度场所需的修正量，同时提升了运行时性能、计算精度和稳定性。
子单元与父单元之间的关系以refinementHistory形式维护，该数据结构同时存储每个单元的加密等级。仅当所有兄弟单元都满足解加密条件时才会执行解加密操作，这一逻辑由dynamicRefineFvMesh类负责处理。否则，在加密区域边界的单元等级可能会在每次加密和解密操作时反复切换，造成不必要的计算资源消耗。
2.2 四叉树细化
对于二维情况，OpenFOAM采用三维网格，并附加了在不变方向上单层单元格的约束条件。在运行时，系统会自动检测到empty边界类型作为二维情况的标识，从而减少需要求解的维度数量。八叉树细化会打破这种不变性条件，因为它在不变方向上引入了新的单元格方向。对于二维情况，标准八叉树细化必须被四叉树细化所取代。四叉树细化的实现基于八叉树细化，不同之处在于内部单元不会沿不变方向分裂产生新单元。空边界上的两个面会在各自中心处分裂为四个新面，这一步骤会在所有边界边的中点引入新顶点，用于沿允许方向分割面。因此单元细化会生成四个子单元，并在父单元内引入四个内部面（见图1(b)）。
对于轴对称情况，该实现采用另一种边界类型——楔形边界。网格仍保持单层单元厚度，但呈现为沿对称轴延伸的楔形结构。如前所述，四叉树细化可应用于六面体单元。但需要特别注意对称轴上棱柱单元的细化过程。这些单元的三角面在细化时会被检测并分割为三角形和四边形面，从而使棱柱单元细分为两个棱柱和两个六面体（见图1(b)）。
将三维网格切割器扩展至二维和轴对称场景的基础来自文献[15]的研究，本研究对其进行了进一步抽象与改进。为减少代码重复，hexRef8类被替换为基础hexRef类及其三个派生类：hexRef8、hexRef4和hexRef4Axi，分别实现八叉树细化、四叉树细化以及包含对称轴上棱柱的四叉树细化（见图2）。hexRef类利用运行时选择机制，使dynamicRefineFvMesh能根据网格维度和求解维度自动选择合适的N叉树细化方案。运行时还会通过维度数判断是否允许进行棱柱单元细化。
2.3 结合多重细化标准
为了获得更大的灵活性，OpenFOAM的单准则细化被扩展为多种不同准则。这些准则包括可单独选择的场、其梯度和旋度、界面等。诸如方盒或域边界等几何特征，以及每个独立准则上的最大与最小细化层级。所有设置均可在运行时修改，其使用方式已在代码库的教程中示范性展示。图2展示了作为动态细化网格模块化插件存在的多准则细化类。
2.4 缓冲层
缓冲层是指两个细化级别之间的单元格层数。在具有两个以上细化级别的网格中，确保不同级别之间的平滑过渡至关重要，这能有效减少因网格在细化过渡处的畸变导致的离散化误差，并为计算流场提供适应新网格级别的缓冲空间。根据”单级不规则约束”规定，相邻单元格的细化级别差异不得超过一级，因此最小缓冲层数为一层。但实际应用中建议设置更多缓冲层。该功能的具体实现方式详见附录B。
3 动态负载均衡
在OpenFOAM中，并行计算基于区域分解法实现。计算域被分解为若干子区域，每个子区域分配给不同的处理器。当在并行瞬态问题中使用自适应网格加密（AMR）时，处理器负载的动态变化可能显著降低计算资源的使用效率。负载最大的处理器会成为性能瓶颈，因为其他等待中的处理器需要其处理完成才能继续工作。
图3 网格细化后的负载均衡再分配。
负载均衡的评估可采用多种方法。本文采用的计算方式是：取各子区域网格数与全局平均网格数的最大差值，再除以该平均值。这种方法假设每个网格单元的计算成本相同，但实际情况往往并非如此。图3展示了一个负载不均衡的案例：当对某组单元进行加密后，区域分布出现56%的不均衡度，这个数值显然过高；经过重新分配后，不均衡度降至仅4%。
目前已有多个基于OpenFOAM库的动态负载均衡（DLB）算法投入使用，例如适用于2.1.x版本的[21]中所述算法，以及适用于2.3.x版本、基于Voskuilen解决方案的中算法—后者也是本研究的算法基础。这些算法都依托OpenFOAM的网格重分配功能，该功能能在给定新旧区域分解方案的情况下，在不同处理器间交换网格单元及其场量数据。如Voskuilen所述，该算法仅需微小改动即可支持带物理量纲场量的重分配。我们在此基础上进一步扩展，增加了表面场量翻转的修正方案，并实现了与二维加密相结合的附加功能。
除Voskuilen列出的问题外，我们在fvMeshDistribute模块中还发现了与第2.1节所述通量翻转问题相关的缺陷—当在不同处理器间重分配网格时，非通量表面场量会被错误翻转，导致求解器崩溃。此外，AMR网格更新时默认的边界值映射方案在与DLB结合使用时存在不足。附录C给出了针对具体应用的解决方案。
3.1 域分解策略
OpenFOAM中支持结构化和非结构化区域分解策略。典型示例包括简单的结构化分解方法（将网格按用户定义的子域数量沿各方向细分）以及非结构化的scotch（并行环境下为ptscotch）分解法（旨在最小化处理器间边界尺寸，参见OpenFOAM用户指南[23]和[24]）。scotch（ptscotch）分解在重新分配时倾向于从头生成新域，这会导致处理器间大量单元数据交换和较高的内存消耗。
无论采用何种分解策略，细化后的兄弟单元不应被分配至不同子域以保持逆细化能力。通过引入新型集群分解方法突破了该限制。自OpenFOAM 4.x版本起，分解约束允许在用户指定的decomposeParDict字典中选择refinementHistory类，从而强制单元族保留在同一子域，使得在保持逆细化能力的同时可采用任意分解策略。
图4 自适应网格细化、动态负载平衡和通量校正（若当前模拟时间与细化间隔匹配）
图4总结了自适应网格细化（AMR）、动态负载平衡（DLB）及通量校正的流程。代码库中的单元测试案例通过最多16个单元的简单配置，对比了AMR和DLB过程中通量翻转修改前后以及表面场增强映射算法的效果差异。
4 结果
本文通过两个典型案例展示AMR（自适应网格加密）和DLB（动态负载平衡）技术在二维与三维模拟中带来的性能提升及计算资源优化效果。三维案例为带障碍物的溃坝模拟，二维案例则为毛细管上升现象。所有计算均采用主频2.5GHz的英特尔®至强®E5-2680 v3处理器，通过interDyMFoam求解器实现，该求解器采用代数型流体体积法捕捉界面。
图5 界面处采用两级细化的溃坝模拟结果，显示t=0.24秒时刻与初始状态
该方法的核心优势在于：在保持处理器间计算负载动态再平衡所实现的并行效率前提下，显著减少了所需网格单元数量。溃坝流动模拟研究不仅有助于理解灾难性溃坝事件机理、提升大坝安全性，还可应用于该案例作为多个计算流体力学(CFD)求解器的验证基准。计算域为边长1米的立方体，中心放置0.4米×0.4米×0.25米的障碍物。如图5所示，水体初始化为0.6米×0.19米×0.75米的长方体。研究对比了两套网格：一套采用空间分辨率3.06·10⁻³米的均匀静态网格，共生成3360万网格单元；另一套在界面区域采用两级自适应网格加密(AMR)，分辨率与静态网格保持一致，总网格量仅260万单元。模拟过程分别启用和禁用动态负载平衡(DLB)功能，允许的最大负载不平衡率设为20%。由于经过同等解析的界面主导了流场特征，自适应网格与非自适应网格的模拟结果未见显著差异。
三种不同运行情况下的模拟加速效果如图6所示。采用自适应网格加密（AMR）和动态负载均衡（DLB）技术后，在使用8、16及125个处理器时均实现了一个数量级的加速。但由于处理器数量增加带来的通信开销，仅获得了次线性加速效果。需注意的是，对于采用AMR且网格单元较少的模拟设置，其最佳处理器数量会相应减少。
为阐明AMR应用中DLB的必要性，研究对比了启用与禁用DLB的模拟情况。使用125个处理器时，采用DLB的模拟耗时较未平衡模拟缩短了一半。可以论证的是，在高度动态的模拟场景中（如关注区域在整个计算域内迁移），由负载不平衡引发的瓶颈效应可能导致并行性能甚至低于静态网格加密方案。
附录D展示了模拟平板间毛细上升的二维案例。与溃坝案例类似，该案例不仅展现出显著加速效果，还证明了二维情况下四叉树加密相较八叉树结构的优势。附录E则对比了汽车空气动力学工业案例中，单独使用AMR与结合负载均衡（LB）的AMR模拟结果。
代码库中附带了上述两个案例的教程，以及管轴对称毛细上升案例。这些教学案例仅限于interDyMFoam求解器的应用，但需说明的是，AMR和DLB功能可通过有限体积法（FVM）应用于OpenFOAM绝大多数流动求解器。
6 结论和前景
我们提出了一种面向对象的方法，用于在OpenFOAM中实现负载均衡的二维和三维自适应网格细化。本文阐述了软件设计、代码结构以及针对工程应用软件可用性的改进方案。我们报告了为提升稳定性和运行时性能所必需的若干库改进。
具体而言，针对OpenFOAM 4.x和5.x版本中六面体拓扑网格，我们设计了二维自适应网格细化方案——该软件此前仅支持三维细化。我们的二维网格细化包含经典二维网格和轴对称网格的细化规则，这些规则通过进一步抽象化的OpenFOAM网格切割器类实现。该开发基于文献[15]中Baniabedalruhman的代码，我们对其进行了增强并抽象为网格切割器类以避免代码重复。我们还开发了模块化框架以支持多种细化准则。
通过改进Voskuilen的动态负载均衡代码，并将其与我们开发的自适应网格细化方案相结合，实现了显著的加速效果。研究表明，仅当结合动态负载均衡时，通过自适应网格细化减少总体单元数量所带来的性能提升才能得以保持。这些改进成果可直接应用于OpenFOAM瞬态顶层求解器处理不可压缩流动问题，其模块化软件设计使得这些功能可以立即投入使用。中短期内，开发工作将致力于实现更高层次的抽象—将动态负载均衡功能从自适应网格细化的具体应用场景中分离出来，从而在模拟过程中实现与OpenFOAM现有动态网格类的协同工作，支持负载均衡的域分解。