PETSc 真的能让 OpenFOAM 如虎添翼吗？

PETSc 是美国阿贡国家实验室开发的轻便式可拓展科学计算工具箱（Portable Extensible Toolkit for Scientific computing），读作 pet-see。将 PETSc 植入 OpenFOAM 替换其原生线性求解器的研究早已有之，但实际效果比较有争议，有研究认为 PETSc 远胜于 OpenFOAM 的原生求解器，也有研究认为 PETSc 在很多情况下远逊于 OpenFOAM 的原生求解器，甚至会花费数倍时间。

这些争议引起了小编的关注，小编决定把玩一下这个大名鼎鼎的工具箱。OpenFOAM® v2006 将 PETSc 作为外部求解器整合进了 ThirdParty 包^[1]，写好了现成的编译脚本，还带有教程。

OpenFOAM® v2006 的编译过程非常容易，如果要使用 PETSc，可以先编译 PETSc，再编译 OpenFOAM。如果你已经有了编译好的 OpenFOAM® v2006，也不用清除，编译好 PETSc，再运行 OpenFOAM 的 Allwmake 脚本就可以更新相应的组件。编译 PETSc 时，容易遇到的问题是网络不好，导致一些源码下载失败。这种情况可以开代理服务器，手动下载，并相应地修改脚本。PETSc 的编译过程参见[1]的最后一节，OpenFOAM 的编译过程参见 OpenFOAM® Quick Build Guide^[2]。

编译好 PETSc 后，运行 PETSc 的示例，位置是 OpenFOAM-v2006/modules/external-solver/tutorials，涉及到 laplacianFoam 和 simpleFoam 两个求解器，但 PETSc 是通用的，也可以在其他求解器中调用。

我们重点研究下 simpleFoam/pitzDaily 这个算例，无脑跑一下 PETSc，再和原生的对比，计算结果非常惊人。PETSc 只用了 1 s，而原生的 pitzDaily 居然用了 10 s。 PETSc 的 log 文件大概是这样的

Time = 1

Initializing PETSc
PETSc-bicg:  Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 1
PETSc-bicg:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 1
PETSc-cg:  Solving for p, Initial residual = 0, Final residual = 0, No Iterations 154
time step continuity errors : sum local = 1.71696, global = -0.010057, cumulative = -0.010057
PETSc-bicg:  Solving for epsilon, Initial residual = 0, Final residual = 0, No Iterations 2
bounding epsilon, min: -1.45735 max: 1080.25 average: 48.855
PETSc-bicg:  Solving for k, Initial residual = 0, Final residual = 0, No Iterations 3
ExecutionTime = 0.57 s  ClockTime = 1 s
Code language: PHP (php)

自带的算例毕竟是海专精的结果，如果把类似的配置无脑迁移到其他算例，是否还能有如此惊人的效果呢？

接下来，我们看一下两种求解器的配置参数，然后迁移到 compressibleInterFoam/depthCharge3D 这个算例。

配置有 PETSc 的 simpleFoam 的 fvSolution 中的 solvers 字段如下

solvers
{
    p
    {
        solver petsc;

        petsc
        {
            options
            {
                ksp_type  cg;
                ksp_cg_single_reduction  true;
                ksp_norm_type none;

// With or without hypre
#if 0
                pc_type  hypre;
                pc_hypre_type boomeramg;
                pc_hypre_boomeramg_max_iter "1";
                pc_hypre_boomeramg_strong_threshold "0.25";
                pc_hypre_boomeramg_grid_sweeps_up "1";
                pc_hypre_boomeramg_grid_sweeps_down "1";
                pc_hypre_boomeramg_agg_nl "1";
                pc_hypre_boomeramg_agg_num_paths "1";
                pc_hypre_boomeramg_max_levels "25";
                pc_hypre_boomeramg_coarsen_type HMIS;
                pc_hypre_boomeramg_interp_type ext+i;
                pc_hypre_boomeramg_P_max "1";
                pc_hypre_boomeramg_truncfactor "0.3";
#else
                pc_type bjacobi;
                sub_pc_type ilu;
#endif
            }

            caching
            {
                matrix
                {
                    update always;
                }

                preconditioner
                {
                    update always;
                }
            }
        }

        tolerance       1e-06;
        relTol          0.1;
    }

    "(U|k|epsilon|omega|f|v2)"
    {
        solver petsc;

        petsc
        {
            options
            {
                ksp_type bicg;
                pc_type bjacobi;
                sub_pc_type ilu;
            }

            caching
            {
                matrix
                {
                    update always;
                }

                preconditioner
                {
                    update always;
                }
            }
        }

        tolerance       1e-05;
        relTol          0.1;
    }
}Code language: PHP (php)

而 OpenFOAM 原生的 fvSolution 中的 solvers 字段如下

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
    }

    "(U|k|epsilon|omega|f|v2)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }
}Code language: JavaScript (javascript)

可见，PETSc 使用 cg 方法、OpenFOAM 默认使用 GAMG 方法来求解压强场。PETSc 提供了非常丰富的求解控制参数。字典中提到的 hypre 是另一个适用于大规模并行的求解器，常常和 PETSc 配合使用，贫民窟男孩并没有大规模并行的实验条件，所以没有编译 hypre。

我将上述配置迁移到 compressibleInterFoam/depthCharge3D，统一了 tolerance。结果真的反转了，PETSc 用的时间更多，用了 7325 s，原生 OpenFOAM 用了 5294 s。可能是我的使用经验不足，没能发挥出 PETSc 的性能，我再调调参数，有了新进展向大家报告。

欢迎有经验的朋友在评论区聊聊心得，也可加微信拍砖（微信号：goodbountyhunter）；不推荐新手入坑 PETSc，这并不是一个可以无脑硬上的算法。

参考

1. ^OpenFOAM® v2006: New and improved numerics https://www.openfoam.com/releases/openfoam-v2006/numerics.php
2. ^OpenFOAM® Quick Build Guide https://develop.openfoam.com/Development/openfoam/blob/develop/doc/Build.md