标准语言并行化：用do concurrent实现海洋模型CPU/GPU统一加速

1. 项目概述：当海洋模型遇上标准并行化

在海洋预报和灾害预警领域，时间就是生命。一场风暴潮的模拟，如果因为计算缓慢而延误数小时，其后果可能是灾难性的。传统的海洋数值模型，如我们团队开发的A2D二维非结构网格模型，虽然物理机制完备，但在面对动辄数万网格、数十万时间步的精细化预报任务时，其串行版本的计算耗时常常成为业务化应用的瓶颈。过去，为了榨取硬件性能，我们不得不在MPI、OpenMP、CUDA甚至OpenACC之间做出选择，每种工具都意味着一套独特的语法、一套额外的学习成本和一份独立的代码维护分支。更棘手的是，为CPU优化的代码往往无法直接跑在GPU上，反之亦然，这导致了巨大的开发和维护开销。

直到我们遇到了标准语言并行化。这并非一个全新的工具，而是现代Fortran语言标准的一部分，核心是一个看似简单的语法结构：do concurrent。它的魅力在于，你写的代码是标准的、符合规范的Fortran，但编译器（如Intel的ifort/ifx或NVIDIA的nvfortran）能读懂你的并行意图，并自动将其转化为针对多核CPU或NVIDIA GPU的高效并行指令。这意味着，同一份源代码，无需任何平台特定的指令，仅通过切换编译器选项，就能在CPU集群或GPU加速器上运行。这听起来像是个美好的愿景，而我们这次的工作，就是将它变成现实，应用于A2D海洋模型，并亲眼验证其效能。

2. 核心思路与技术选型：为什么是do concurrent？

在深入代码之前，我们得先理清一个根本问题：在众多并行工具中，为何独独青睐do concurrent？这背后是对开发效率、维护成本和长期可移植性的综合考量。

2.1 传统并行方案的困境

让我们快速回顾一下主流并行编程模型的优缺点：

• MPI (Message Passing Interface)：分布式内存并行的事实标准，擅长跨节点、跨机器的超大规模计算。但其“消息传递”的范式要求对计算域进行复杂的网格分区（Domain Decomposition），代码侵入性强，调试困难，且主要解决的是“多机”并行，对单节点内“多核”或“众核”的利用并非其专长。
• OpenMP：共享内存并行模型，通过编译指导语句（如!$omp parallel do）实现多核CPU并行。语法相对简单，但仅限于CPU，无法利用GPU的庞大算力。在当今GPU加速计算普及的背景下，其能力显得局限。
• CUDA：NVIDIA GPU编程的利器，性能天花板最高。但代价是代码完全与NVIDIA硬件绑定，需要显式管理设备内存、线程网格、块等概念，学习曲线陡峭，代码可读性和可移植性极差。
• OpenACC：通过编译指导语句实现GPU加速，比CUDA更友好，但依然是指令式的，代码中遍布!$acc开头的注释行。虽然支持多平台，但不同编译器的支持程度和细节常有差异。

核心矛盾在于：为了兼顾CPU和GPU，你往往需要维护两套甚至三套代码逻辑相似但语法迥异的源码分支。任何算法修改都需要在多处同步，极易出错，维护成本呈指数级增长。

2.2do concurrent的优势与原理

do concurrent的出现，正是为了化解这一矛盾。它不是一个外部库，而是Fortran 2008及后续标准中定义的语言特性。其基本形式如下：

do concurrent (i = 1:n) a(i) = b(i) + c(i) * d(i) end do

你告诉编译器：这个循环的每次迭代是独立的，没有数据依赖（即计算a(i)时不会读取或写入a(j)，其中j≠i）。基于这个承诺，编译器可以自由地以任何顺序、在任何可用的并行硬件上执行这些迭代。

它的技术价值体现在：

1. 代码简洁与可读性：代码中没有任何!$omp或!$acc指令，纯粹是标准Fortran。任何熟悉Fortran的开发者都能轻松阅读和维护。
2. 硬件抽象与可移植性：并行化的责任从开发者转移到了编译器。开发者只需关心“什么可以并行”，而“如何并行”则由编译器根据目标硬件（多核CPU或GPU）和编译选项来决定。
3. 降低维护成本：单一代码库支持多种硬件目标。无论是为实验室的GPU工作站优化，还是为超算中心的CPU集群部署，都无需修改核心计算代码。
4. 与现代化工具链集成：主流HPC编译器（Intel oneAPI的ifort/ifx， NVIDIA HPC SDK的nvfortran， GNU的gfortran也在逐步完善支持）都积极支持该特性，并利用其底层的OpenMP或OpenACC运行时来实现加速，享受成熟生态的支持。

当然，天下没有免费的午餐。do concurrent将性能优化的部分控制权交给了编译器，其生成的代码性能可能不如手工精心调优的CUDA内核。但对于像A2D模型这样包含大量规整循环（遍历网格单元、边、节点）的科学计算程序，编译器已经能够生成非常高效的代码，在开发效率与运行效率之间取得了极佳的平衡。

3. A2D模型并行化改造实战

理论很美好，实践是检验真理的唯一标准。我们的战场是A2D模型的源代码。这是一套用Fortran编写的、基于有限体积法的二维海洋动力学模型，用于模拟风暴潮、潮汐等过程。

3.1 模型结构与并行潜力分析

首先，我们对模型进行了“体检”，识别计算热点。模型的主循环是一个时间推进过程，在每个时间步内，依次调用一系列子程序来计算水位、流速、通量等。通过分析，我们发现超过95%的计算时间都消耗在几个核心的子程序中，而这些子程序的核心又是对网格数组（如所有单元、所有边）进行遍历的do循环。

例如，计算水平通量的一个典型循环原先是这样的：

do j=1, jm hmflux(j) = ua(j) * ds(j) * xn(j) + va(j) * ds(j) * yn(j) end do

这里，jm是边的数量，hmflux、ua、va、ds、xn、yn都是一维数组。这个循环对每个边j进行独立计算，完美符合do concurrent“无数据依赖”的前提。

这里有一个至关重要的细节：A2D模型虽然使用二维三角形非结构网格，但在内存中，所有网格相关的变量（单元中心的水位、边中心的流速等）都存储在一维数组中。这种“一维化”的存储方式，使得我们的并行化改造变得异常简单和规整，只需要处理一维循环，避免了多维循环可能带来的循环次序优化等复杂问题。

3.2 并行化改造的具体步骤

改造的核心工作，就是将模型中所有可并行的do循环替换为do concurrent。以上述通量计算为例，改造后的代码如下：

do concurrent (j=1:jm) local(xmflux, ymflux) xmflux = ua(j) * ds(j) ymflux = va(j) * ds(j) hmflux(j) = xmflux * xn(j) + ymflux * yn(j) end do

关键改动与解释：

1. 循环头替换：do j=1, jm->do concurrent (j=1:jm)。这是最主要的语法变化。
2. local子句：我们显式声明了xmflux和ymflux这两个标量变量为local。这意味着每个并行迭代（例如每个GPU线程或CPU线程）都会拥有自己独立的这两个变量的副本，从而避免多个迭代同时读写同一个内存地址导致的“竞态条件”。虽然一些编译器在简单情况下能自动推断，但显式声明local是保证代码在不同编译器间行为一致的最佳实践，也是写出安全并行代码的好习惯。
3. 数据依赖检查：这是改造过程中最需要谨慎的一步。我们必须确保循环体内，对数组的写操作（如hmflux(j) = ...）的索引j，与右侧任何数组的读索��无关。例如，绝不能出现hmflux(j) = hmflux(j-1) + ...这样的情况，因为这会引入迭代间的依赖，破坏并行安全性。A2D模型的算法本身是显式且局部性很好的，这为我们省去了大量重构算法的工作。

实操心得一：改造不是简单的“查找替换”

• 全局搜索与人工审核：我们使用文本工具全局搜索所有的do循环，但并非所有循环都适合并行。需要仔细审查循环体，识别真正的数据依赖。例如，包含sum = sum + a(i)这类归约操作的循环，需要特殊处理（标准Fortran 2018后do concurrent支持reduce子句，但需要编译器支持）。
• I/O与边界条件：模型中确实存在无法并行化的部分，主要是文件输出(OutputVd等子程序)和依赖于时间序列的边界条件更新(BCOND_Time中的部分代码)。这些是固有的串行部分，我们保留其原貌。根据阿姆达尔定律，这些部分将决定并行加速的上限。
• 数组分配方式：一个容易被忽略但至关重要的点是，要使代码能运行在GPU上，所有需要在GPU上参与计算的数组必须是动态分配的（allocatable）。这是因为NVIDIA的Unified Memory（统一内存）技术依赖于动态内存管理来实现CPU和GPU之间的数据自动迁移。固定大小的静态数组（dimension(1000)）无法享受此便利。幸运的是，A2D模型原本就大量使用了可分配数组。

3.3 编译与运行：一键切换硬件平台

改造完代码后，最激动人心的时刻到来：通过编译器选项来控制执行平台。

• 为多核CPU编译 (使用Intel编译器):

  ifort -qopenmp -O2 a2d_parallel.f90 -o a2d_cpu.exe

  这里，-qopenmp告诉Intel编译器，将do concurrent通过OpenMP后端实现。运行时可使用环境变量OMP_NUM_THREADS或代码内调用call omp_set_num_threads(N)来设置线程数。

• 为多核CPU编译 (使用NVIDIA编译器):

  nvfortran -stdpar=multicore -O4 a2d_parallel.f90 -o a2d_cpu_nv.exe

  -stdpar=multicore指示nvfortran在CPU上执行标准并行代码。-O4是NVIDIA编译器的高级别优化选项，实测能带来小幅性能提升。

• 为NVIDIA GPU编译:

  nvfortran -stdpar=gpu a2d_parallel.f90 -o a2d_gpu.exe

  -stdpar=gpu是魔法发生的地方。编译器会自动将do concurrent循环卸载到GPU上执行，并利用Unified Memory自动管理CPU和GPU之间的数据传输，开发者无需手动拷贝数据。

实操心得二：编译器的“小脾气”不同编译器对do concurrent的实现和优化策略有细微差别。例如，在循环嵌套的情况下，do concurrent (i=1:im, j=1:jm)和do concurrent (j=1:jm, i=1:im)的性能可能不同，这取决于编译器如何映射迭代到线程块。由于A2D模型采用一维数组存储，我们幸运地避开了这个复杂性。但如果你正在处理多维循环，建议进行简单的性能测试来确定最佳循环次序。

4. 性能测试与结果分析

我们在两个高性能计算平台上进行了测试：

• 平台1：配备双路28核Intel Xeon Gold CPU和两张NVIDIA A800 GPU的节点。
• 平台2：配备双路32核Intel Xeon Platinum CPU的节点。

我们以原始的串行版本为基准，对比了不同核心数下的CPU并行版本以及GPU版本的加速效果。每个配置都运行多次取平均，以减小误差。

4.1 多核CPU性能表现

测试结果令人振奋。在平台1（56核）上，无论是使用ifort还是nvfortran，并行版本都取得了超过23倍的加速比。在平台2（60核）上，加速比更是超过了25倍。下图直观展示了随着核心数增加，计算时间迅速下降的趋势。 （注：此处原应为图表，在Markdown中以文字描述代替）执行时间随CPU核心数增加而近乎线性下降，在56核时降至串行版本的约1/25。

理论极限与实际情况：我们根据精细的性能剖析，测量了代码中固有的串行部分比例（α）。在单核情况下，这个比例极小（约0.05%~0.1%）。根据阿姆达尔定律，即使串行部分只有0.1%，在100个处理器上的理论加速上限也在90倍以上。我们实测的20-30倍加速，说明除了串行部分，还受到了内存带宽、缓存一致性、线程同步开销等因素的限制。但即便如此，20倍以上的性能提升对于实际业务预报而言，意味着将原本需要10小时的模拟缩短到30分钟以内，这是质的飞跃。

4.2 GPU卸载性能与一个“坑”

使用nvfortran -stdpar=gpu编译运行在A800 GPU上，我们最初遇到了一个有趣的现象：使用编译器版本23.7时，获得了超过23倍的加速（约14.69秒 vs 串行344.13秒），但与版本24.7和25.7时，性能却出现了倒退。

问题排查：我们通过添加编译器诊断标志-Minfo=acc来查看编译器背后做了什么。发现新版本的编译器对一些嵌套在do concurrent内部的短循环（例如，循环上限只有8的小循环）自动尝试了“隐式归约”优化。这本是好事，但对于这种迭代次数极少的循环，启动并行带来的开销远大于收益，反而拖累了整体性能。

解决方案：对于这5个特定的、被编译器“过度优化”的循环，我们没有放弃标准并行化。而是采用了一种混合策略：将这5个循环改回传统的do循环，但为其加上明确的OpenACC指令!$acc loop seq，告诉编译器“这个循环请保持串行”。同时，将编译选项改为-stdpar=gpu -acc=gpu。

! 原始 do concurrent 结构（可能被过度优化） do concurrent (m=1:mb) do n=1, ns_n(m) ! 这个内层循环很小 ! ... 计算 ... end do end do ! 修改为混合模式 do concurrent (m=1:mb) !$acc loop seq do n=1, ns_n(m) ! ... 计算 ... end do end do

经过这样微调后，再用nvfortran 24.7/25.7编译，性能恢复到了与23.7相当的水平，加速比稳定在23倍以上。

实操心得三：拥抱标准，但不排斥工具这个插曲说明了两个问题：1) 编译器的优化策略在不断演进，有时会引入非预期的性能变化；2) 标准语言并行化与传统的指令式并行（如OpenACC）并非互斥。在绝大多数循环使用do concurrent保证简洁和可移植性的前提下，对极少数特殊场景使用精准的OpenACC指令进行微调，是一种务实且高效的策略。最终，我们仍然保持了代码主体的高度可读性和可移植性。

4.3 数值一致性验证

并行计算，速度很重要，但正确性更重要。我们对比了并行版本与串行版本在10天模拟期内，关键潮位站（如岱山、吕四）的输出结果。计算了逐日的均方根误差。

结果：所有并行配置（8核、28核、56核CPU，以及GPU）与串行参考解之间的差异微乎其微。在整个10天的模拟中，水位结果的RMSE最大不超过4毫米（岱山站），吕四站甚至小于0.8毫米。考虑到风暴潮模拟中水位变幅通常在半米到两米之间，这样的误差完全在可接受的数值精度范围内，对预报结果没有实质影响。这强有力地证明了我们的并行化改造没有引入算法错误，数值结果是可信的。

5. 经验总结与未来展望

回顾整个项目，从分析、改造、调试到测试，我们验证了标准语言并行化在海��数值模型加速中的可行性与高效性。以下是一些可供同行参考的经验：

1. 从“热点”开始，循序渐进：不要试图一次性并行化所有代码。先用性能分析工具（如gprof,nvprof）找到最耗时的子程序，优先改造其中的循环。成功后再逐步扩大范围。
2. 严格审查数据依赖：这是并行化安全的重中之重。对于每一个do concurrent循环，必须像代码审查一样，确认没有跨迭代的读写冲突。利用编译器的检查功能（如Intel编译器的-qopt-report-phase=loop）可以提供帮助。
3. 善用编译器的诊断信息：无论是Intel的-qopt-report还是NVIDIA的-Minfo=acc，这些输出信息能让你看清编译器背后做了什么优化、哪些循环被并行化了、数据是否被传输到设备等，是调试和性能调优的利器。
4. 理解硬件特性：CPU并行和GPU并行有不同的优化侧重点。CPU核心少，但单核能力强，适合复杂逻辑；GPU核心多，但单核弱，适合大规模数据并行。do concurrent虽然做了抽象，但在编写循环体时，潜意识里考虑目标硬件，避免在循环内进行大量分支判断或调用复杂函数，有助于编译器生成更高效的代码。
5. 混合并行是未来：对于超大规模的模拟，单个计算节点（无论是多核CPU还是多GPU）可能仍不够。未来的方向可能是“MPI + 标准语言并行化”的混合模式。MPI负责跨节点分布式并行，而每个节点内部则用do concurrent充分利用所有CPU核心或GPU。这样既能扩展至成千上万个节点，又能高效利用节点内资源。

最后一点个人体会：标准语言并行化，特别是do concurrent，它降低的不仅仅是代码的维护成本，更是一种思维负担的减轻。开发者可以更专注于科学问题本身和算法的实现，而将如何高效利用并行硬件的问题，更多地交给日益智能的编译器工具链。这对于海洋、大气、地质等领域的广大科研人员和工程师来说，无疑是一大利好。它让高性能计算的门槛变得更低，让更复杂的模型、更高分辨率的模拟、更快速的业务化预报成为可能。我们的A2D模型实践只是一个开始，期待看到这项技术在未来更广泛的领域开花结果。