自定义前处理应用程序可能涉及以 OpenFOAM 中尚不可用的方式准备初始条件、准备在 HPC 集群上执行并行模拟或设置参数变化。在本节中,我们关注并行执行和参数变化。

使用通过 shell 或基于 Python 的脚本链接在一起的各种应用程序,通常可以执行和自动化数量惊人的任务。 OpenFOAM 的良好做法是在模拟案例文件夹中准备一个所谓的 Allrun 脚本,该脚本执行准备模拟所需的预处理脚本。

本节介绍了两个具有不同复杂性的示例。在第一个示例中,一个简单的 shell 脚本使用现有的预处理应用程序调用 OpenFOAM 可执行文件。第二个示例涵盖使用 PyFoam 库对预处理应用程序进行编程以生成参数变化。 Bernhard Gschaider 主要开发 PyFoam 项目:一组用 Python 编程语言编写的库和可执行程序,用于直接参数化和分析 OpenFOAM 模拟。

8.2.1 并行应用程序执行

考虑大量模拟需要使用高性能集群的情况。对于这个例子,所有的模拟案例将被分解成不同的子域,并存储在一个名为 Simulations 的目录中。模拟目录将作为 $FOAM_RUN 目录的子目录放置。与其手动启动每个模拟,不如尽可能自动化此过程。为此,可以使用以下步骤对 shell 脚本进行编程:

  1. 从 controlDict 检索求解器名称
  2. 从 decomposeParDict 中检索子域的数量。
  3. 使用正确数量的处理器执行 mpirun 命令。

在此示例中,脚本名为 Allparrun,并存储在模拟目录中。作为第一步,应该在脚本的开头添加几行。在随 OpenFOAM 发布的大部分教程中都可以找到类似的代码:

#!/bin/bash
cd ${0%/*} || exit 1
source $WM_PROJECT_DIR/bin/tools/RunFunctions
application=getApplication

getApplication 函数是 OpenFOAM 中现有的 bash 函数,非常适合本示例,因为在本例中,应用程序的名称将是求解器的名称。

接下来,必须找到simulation目录中的所有案例子目录。假设这些目录都是正确的 OpenFOAM 模拟案例。因此,OpenFOAM 模拟将在模拟目录的每个子目录中启动。 find 命令用于专门查找可在 $FOAM_RUN/simulations 中找到的目录(参见清单 9)。上述行查找直接(非递归)位于 $FOAM_RUN 目录中的任何目录,以便稍后可以对这些情况中的每一个执行相同的操作。 find 命令可以很容易地通过例如使用通配符过滤某些目录名称,但为了简单起见,在此示例中假设模拟目录仅包含 OpenFOAM 案例。

for d in `find $FOAM_RUN/simulations -type d -maxdepth 1`
do
# This part will be programmed at a later point in this section.
done

下一步是从 system/controlDict 检索子域的数量并将其存储在变量中以供以后使用。为此,在 Allparrun 脚本中引入了一个名为 nProcs 的新函数,该函数返回子域的数量:

function nProcs
{
    n=`grep "numberOfSubdomains" system/decomposeParDict \
    | sed "s/numberOfSubdomains\s//g" \
    | sed "s/;//g"`
    echo $n
}

该函数调用 grep 程序来过滤 system/decomposeParDict 中包含 numberOfSubdomains 的行。在这一行中,“numberOfSubdomains”字符串被删除,该关键字和实际数值之间的任何空格都被删除,并且尾分号被删除。此链接命令的结果存储在 bash 变量 n 中。要返回它,使用 echo 。

使用 mpirun 命令并行运行任何 OpenFOAM 求解器或实用程序。对于此示例,这意味着需要将两个变量作为选项参数传递给 mpirun:求解器名称和处理器数量。完成的 Allparrun 脚本如下所示。

#!/bin/bash
cd ${0%/*} || exit 1
source $WM_PROJECT_DIR/bin/tools/RunFunctions
function nProcs
{
    n=`grep "numberOfSubdomains" system/decomposeParDict \
    | sed "s/numberOfSubdomains\s//g" \
    | sed "s/;//g"`
    echo $n
}
PWD=`pwd`

for d in `find $FOAM_RUN/simulations -type d -maxdepth 1`
do
    # Jump into the current case directory
    cd $d
    n=`nProcs`
    application=`getApplication`
    runApplication decomposePar
    # Depending on the environment of your HPC, this must be
    # adjusted accordingly. Remember to provide a proper
    # machine file, otherwise all jobs will be started on the
    # master node, which is undesirable in any way.
    mpirun -np $n $application -parallel >  log &
    # Jump back
    cd $PWD
done

对于最简单的情况,作业排队和提交系统不可用,您只需向 mpirun 提供一个machine文件。machine文件包含 HPC 集群网络上计算节点的 IP 地址或主机名。有关此主题的更多信息,请参阅您的 HPC 集群管理团队提供的相应用户手册。

提示:请记住,在这种情况下,对 mpirun 的调用会将所有作业作为机器上的本地作业启动。为了在 HPC 集群上启动它并使用其计算节点,其他步骤是必要的,具体取决于特定的集群配置,因此此处省略。

8.2.2 参数变化

参数变化包括对具有不同物理和离散化参数的许多模拟的预处理,这需要自动化。为了在 OpenFOAM 中模拟参数研究,可以使用 shell 脚本、Python 脚本或可用的 OpenFOAM 实用程序来克隆模拟并在基于文本的输入(配置)文件中调整它们的参数。由于 OpenFOAM 是为命令行使用而构建的,因此 OpenFOAM 比其他基于 GUI 的 CFD 平台更适合自动化。

此示例涵盖二维 NACA0012 翼型的参数化。网格生成是使用 blockMesh 应用程序完成的,网格使用 mirrorMesh 应用程序进行镜像。参数研究调查了攻角α对升力和流动空气对翼型施加的阻力的影响。参数研究由 15 次模拟组成,攻角从 α = 0° 到 α = 15°,导致升力 w.r.t 的粗略分布。然后通过在前 15 个收敛后自动创建新模拟,在图中的峰值处细化初始升力分布。如果没有自动化,CFD 项目中的大量参数变化会很快成为一项费力且容易出错的活动。

使用 shell 脚本进行自动化是可能的,但很麻烦,因为有效地操作 OpenFOAM 配置文件需要大量的源代码,以及 shell 脚本如何处理字母数字计算。

Python 脚本语言可以直接替代 shell 脚本,进一步简化了模拟结果的处理和可视化。此外,许多基于 Python 的库已经具有旨在支持数值模拟过程的功能,例如插值方法、求根方法、数据处理和分析库、机器学习等。一个 Python 项目已经实现了 OpenFOAM 模拟的参数变化自动化所需的许多功能:PyFoam。

PyFoam 项目由可以操作 OpenFOAM 数据的不同 Python 模块组成。除了这个库部分,PyFoam 已经提供了许多不同的即用型应用程序来简化不同的任务。

命令行中的 Tab-completion 显示了许多可用 PyFoam 应用程序的列表;在这里,我们只评论其中的几个:pyFoamClearCase.py 清除 OpenFOAM 模拟案例目录中先前计算的结果,pyFoamPlotWatcher.py 从 OpenFOAM 求解器的日志文件中实时呈现残差图,pyFoamRunParameterVariation.py 创建参数研究来自参数向量数组和模板案例的笛卡尔积。

PyFoam 可执行文件建立在 PyFoam 库之上,与以下涵盖 PyFoam 预处理应用程序编程的示例相同。在 Python 脚本中使用 PyFoam 库可以对 OpenFOAM 案例进行自定义操作。此外,可以使用 subprocess.call Python 子模块以直接的方式在 Python 中执行 OpenFOAM 可执行文件。对于这个例子,我们假设有一些 Python 的基本背景知识,并且机器上可以使用 PyFoam。

提示:PyFoam 可以使用 Python 的 pip 包安装程序通过 sudo pip install PyFoam 安装,或者在没有 root 权限的情况下使用 pip install --user PyFoam 安装。

尽管可以使用 PyFoam 从 OpenFOAM 操作 CFD 结果,但必须了解其操作大型数据数组(场变量值)的运行速度比使用 C++ 等编译语言要慢一些。这个例子侧重于参数变化研究的准备,而不是处理 OpenFOAM 模拟产生的物理场。

OpenFOAM 中的参数研究由许多具有不同配置文件的模拟文件夹组成。模拟文件夹都是从相同的输入(模板)模拟文件夹生成的。在这个例子中,我们将重新使用第 4 章中的模拟:ofprimer/cases/chapter04/naca 作为参数变化的模板案例。

提示:正确设置模板模拟案例至关重要,因为所有参数化案例都是从中创建的。

在为更复杂的参数研究编写预处理应用程序之前,有必要了解使用 pyFoamRunParameterVariation.py 已经可以实现的目标。 PyFoam 应用程序具有大量文档,可以使用命令行中的 --help 选项显示这些文档。执行 pyFoamRunParameterVariation.py --help 会显示一个广泛的选项列表,此处并未全部涵盖。

将参数研究目录结构的生成与真实场景中的网格划分和仿真执行分开是一种很好的做法。这样,网格生成和模拟在分解问题的子域和研究参数中并行运行。换句话说,我们首先串行生成模拟文件夹结构,因为这仅涉及文本文件的操作。然后使用 HPC 资源上的工作负载管理器对每个参数向量并行运行网格化,然后进行各自的模拟。为了实现这种设置,pyFoamRunParameterVariation.py 需要一组特定的选项,这些选项包含在名为 create_naca_study.py 的 Python 脚本中,以便在有效的工作流程中重复使用。

提示:脚本 create_naca_study.py 可以在 naca 模板模拟案例旁边的 ofprimer/cases/chapter04 文件夹中找到。

在 create_naca_study.py 中传递给 pyFoamRunParameterVariation.py 的选项如下所示:

  • --every-variant-one-case-execution 确保为变量中的每个参数向量生成一个 OpenFOAM 模拟案例,
  • --create-database 存储模拟 ID 和参数向量之间的关系(对于再现结果至关重要)
  • --no-mesh-create 防止网格生成
  • --no-execute-solver 阻止求解器(模拟)执行
  • --cloned-case-prefix 用于将此参数研究与其他研究区分开来

在 create_naca_study.py 脚本(如下所示)中存储对 pyFoamRunParameterVariation.py 的此类调用的“配置”很有用:这简化了模拟工作流程并确保轻松再现结果。

#!/usr/bin/env python
import argparse
import sys
from subprocess import call
parser = argparse.ArgumentParser(description='Runs \
                    pyFoamRunParameterVariation.py to create simulation folders \
                    without generating the mesh.')
parser.add_argument('--study_name', dest="study_name", \
                    type=str, required=True, \
                    help='Name of the parameter study.')
parser.add_argument('--parameter_file', dest="parameter_file", \
                    type=str, required=True,
                    help='PyFoam .parameter file')
parser.add_argument('--template_case', dest="template_case", \
                    type=str, required=True,
                    help='Parameter study template case.')
args = parser.parse_args()

if __name__ == '__main__':
    args = parser.parse_args(sys.argv[1:])
    prefix = args.study_name + "_" + args.parameter_file
    call_args = ["pyFoamRunParameterVariation.py",
            "--every-variant-one-case-execution",
            "--create-database",
            "--no-mesh-create",
            "--no-server-process",
            "--no-execute-solver",
            "--no-case-setup",
            "--cloned-case-prefix=%s" % prefix,
            args.template_case,
            args.parameter_file]
    call(call_args)

下列程序片段显示了 naca 案例的参数,其中只有运动粘度发生变化。其列出并存储在 naca.parameter 文件中的值应该在 naca 案例的 constant/transportProperties 文件中以某种方式替换。

values
{
    solver ( simpleFoam );
    NU
    (
        1e-06 2e-06 3e-06 4e-06
    );
}

为确保 PyFoam 替换 constant/transportProperties 中的 值,将该文件复制到 constant/transportProperties.parameter 中,并修改 nu,如下所示。

...
transportModel Newtonian;
 rho        rho [ 1 -3 0 0 0 0 0 ] 1000;
 nu        nu    [0 2 -1 0 0 0 0] @!NU!@;
...

运行以下命令:

?> ./create_naca_study.py --study_name my_study --template_case naca --parameter_file naca.parameter

从naca模拟案例模板生成清单13中列出的四个具有四个粘度值的OpenFOAM新模拟案例:

?> ls -d my_study*
my_study_naca.parameter_00000_naca
my_study_naca.parameter_00001_naca
my_study_naca.parameter_00002_naca
my_study_naca.parameter_00003_naca

如果清单13中定义了多个参数向量,则PyFoam将构造这些向量及其各自的模拟文件夹的笛卡尔积。大的参数变化要求唯一的模拟id (00000、00001、……) 和参数向量之间的关系。这很容易用下面的命令实现:

?> pyFoamRunParameterVariation.py --list-variations naca naca.parameter
==============================
4 variations with 1 parameters
==============================
==================
Listing variations
==================
Variation 0 : {'NU': 1e-06}
Variation 1 : {'NU': 2e-06}
Variation 2 : {'NU': 3e-06}
Variation 3 : {'NU': 4e-06}

用PyFoam的pyFoamRunParameterVariation.py涵盖了现有的OpenFOAM仿真参数的变化,接下来的步骤涉及网格的生成和仿真。

下一个示例涵盖使用PyFoam库中的子模块对参数变化脚本进行编程。当参数变化的生成超出参数向量的笛卡尔乘积,或者更复杂的计算确定参数值时,此任务变得相关。

在详细介绍编写自定义的PyFoam脚本之前,将提供有关仿真的必要背景信息。使用块网格生成对称网格,该块网格由图8.1中概述的三个块组成。调整块的分级,以使从块1到2和从2到3的网格大小没有明显变化。为了解决沿翼型的粘性边界层,在NACA轮廓附近生成了相当小的单元。

根据图8.1中所示的标记来命名所有patch。顶点数由两个用破折号分隔的数字表示。机翼上未沿图形法线方向定向的所有边均使用blockMesh的polyLine生成,其中根据定义任何两位数NACA轮廓的函数计算站。蓝色入口的顶点也是如此-但不是使用两位数的NACA多项式,而是使用四分之一圆。这使网格看起来更具视觉吸引力,并消除了入口和侧面之间的任何边界条件问题。

图8.1 NACA网格概述

计算域的所有维度以及NACA分布本身都可以从图8.1中获得。 所采用的求解器是simpleFoam,并将以串行方式执行。

注:cases/chapter 04/naca案例仅用于网格生成和参数变化,而不是用于获得物理结果。

Python脚本的编程从导入必要的模块开始,如清单14所示。

# 导入PyFoam包
from PyFoam.RunDictionary.SolutionDirectory import SolutionDirectory
from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
from PyFoam.Basics.DataStructures import Vector
from numpy import linspace,sin,cos,array
from os import path, getcwd

对于本节中介绍的示例,有三个PyFoam子模块非常重要:

  • SolutionDirectory:处理整个案例并提供对其数据的轻松访问以及将案例克隆到不同目录的方法
  • ParsedParameterFile:读取任何OpenFoam字典并提供对它的访问,就像它是Python字典一样
  • Vector:处理向量的操作

用于解析命令行参数和文件路径以及在自定义参数变体中执行计算的模块(参见清单15)直接跟随清单14中的PyFoam模块。

# 程序清单15
# 导入PyFoam模块
# Commandline arguments and paths
import os
import sys
import argparse
# Parameter variation
import numpy as np
from numpy import linspace,sin,cos,array,pi
import yaml

自定义参数研究将改变NACA0012翼型的迎角,在本例中,这需要向量绕-y坐标轴旋转,由清单16中的函数实现。

# 程序清单16
def rot_matrix(angle_deg):
    angle_rad = angle_deg * pi / 180.
    return np.array([[cos(angle_rad),0,-sin(angle_rad)],[0,1,0]
                     [sin(angle_rad), 0, cos(angle_rad)]])

def rot_vector(vec, angle_deg):
    return np.dot(rot_matrix(angle_deg), vec)

参数变化脚本通过命令行选项控制学习。 或者,可以像pyFoamRunParameterVariation.py一样控制参数变化:使用配置文件,采用CSV、JSON或YAML等标准化格式。 命令行参数由清单17所示的代码解析。 参数变化脚本的主要功能如清单18所示,并用注释对其进行了详细概述。 其核心思想是读取一个模板仿真案例,然后将该案例克隆为多个仿真案例,这些案例由与参数变化相对应的唯一标识符标识。

参数变化脚本通过命令行选项控制研究。或者可以pyFoamRunParameterVariation.py使用配置文件(采用CSV、JSON或YAML等标准化格式)来控制参数变化。命令行参数由清单17所示的代码解析。参数变化脚本的主要功能如清单18所示,并通过注释详细概述了它。核心思想是读取模板仿真案例,然后将该案例克隆成多个仿真案例,由对应于参数变化的唯一标识符(ID)来标识。

# 程序清单17
parser = argparse.ArgumentParser(description='Generates simulation cases for parameter study using PyFoam.')
parser.add_argument('--template_case',dest="template_case", type=str,help='OpenFOAM template case.', required=True)
parser.add_argument('--study_name', dest="study_name", type=str,help='Name of the parameter study.', required=True)
parser.add_argument('--alpha_min', dest="alpha_min", type=float,help='Minimal angle of attack in degrees.',required=True)
parser.add_argument('--alpha_max', dest="alpha_max", type=float,help='Maximal angle of attack in degrees.',required=True)
parser.add_argument('--n_alphas', dest="n_alphas", type=int,help='Number of angles between alpha_min and alpha_max.',required=True)
args = parser.parse_args()
# 程序清单18
if __name__ == '__main__':
    args = parser.parse_args(sys.argv[1:])
    # Distribute angles of attack
    angles = np.linspace(args.alpha_min, args.alpha_max, args.n_alphas)
    # Initialize the template case
    template_case = SolutionDirectory(args.template_case)
    # Initialize the parameter dictionary
    param_dict = {}
    # For every variation
    for variation_id,alpha in enumerate(angles):
    # Add variation-d : parameter vector to parameter dictionary.
    param_dict[variation_id] = {"ALPHA" : float(alpha)}
    # Clone the template case into a variation ID case.
    name = "%s-%08d" % (args.study_name, variation_id)
    case_name = os.path.join(os.path.curdir, name)
    cloned_case = template_case.cloneCase(case_name)
    # Read the velocity field file.
    u_file_path = os.path.join(cloned_case.name,"0", "U")
    u_file = ParsedParameterFile(u_file_path)
    # Read the velocity inlet boundary condition and internal field.
    u_inlet = u_file["boundaryField"]["INLET"]["value"]
    u_internal = u_file["internalField"]
    # Rotate the velocity clockwise by alpha around (-y).
    u_numpy = np.array(u_inlet[1])
    u_numpy = rot_vector(u_numpy, alpha)
    # Set inlet and internal velocity.
    u_inlet.setUniform(Vector(u_numpy[0],u_numpy[1],u_numpy[2]))
    u_internal.setUniform(Vector(u_numpy[0],u_numpy[1],u_numpy[2]))
    u_file.writeFile()
    # Store variation ID : study parameters into a YAML file.
    with open(args.study_name + '.yml', 'w') as param_file:
    yaml.dump(param_dict, param_file, default_flow_style=False)

此示例中的每个参数变量都由单个标量值组成:迎角α。通常,通过参数空间中的n维参数向量来识别参数变化,其中n对应于仿真中变化的参数的数目。在这种情况下,改变迎角需要速度场围绕−y坐标轴方向旋转。

使用PyFoam读取速度文件,并从中获得入口边界条件和内部速度场。然后设置入口边界条件和初始内部值,以达到特定的迎角。一旦像这样初始化速度场,就会写入其文件。然后将唯一参数变量的ID与参数向量相关联的字典保存为YAML文件。这个ID -参数向量映射也由pyFoamRunParameterVariation.py以PyFoam特定的格式保存。这里,YAML被用作一种标准格式,在Python中很容易操作。尽管在这个小的示例脚本中只有一个参数发生了变化,但该示例显示了如何将PyFoam的主要子模块与标准Python模块一起使用,以扩展OpenFOAM中的参数变化,从而完全控制该变化。例如,下一步是通过应用条件来禁用没有意义的参数向量,从而减少参数向量的全笛卡尔积。

该脚本位于代码库中,如cases/chapter 04/parametricStudy.py。称之为在0 °和10 °之间以10个步长改变迎角,即

?>  ./parametricStudy.py --study_name angle-of-attack \
--template_case naca --alpha_min 0 \
--alpha_max 10 --n_alpha 10

得到结果为:

?>  ls
angle-of-attack-00000001
angle-of-attack-00000002
angle-of-attack-00000003
...
angle-of-attack.yml

angle-of-attack.yml文件中包含ID参数向量映射:

0:
ALPHA: 0.0
1:
ALPHA: 1.1111111111111112
2:
ALPHA: 2.2222222222222223
...

虽然naca情况的模拟参数不能保证本例中物理结果的正确性,但参数研究的两个极端变量的流线如图8. 2所示。

图8.2 α = 0(迎角= 00000000)和α = 10(迎角= 00000009)时的流线