# palabosUserGuide **Repository Path**: zhangqingdian/palabosUserGuide ## Basic Information - **Project Name**: palabosUserGuide - **Description**: 英文版palabos用户手册翻译 - **Primary Language**: Unknown - **License**: Not specified - **Default Branch**: master - **Homepage**: None - **GVP Project**: No ## Statistics - **Stars**: 2 - **Forks**: 1 - **Created**: 2021-12-27 - **Last Updated**: 2024-03-04 ## Categories & Tags **Categories**: Uncategorized **Tags**: None ## README # Palabos用户手册中文版 本用户手册翻译自Palabos User Guide(Release 1.0 r1)。 译者水平不高，大部分翻译来自机翻，会有很多语句不通顺的地方；另外，由于译者所涉猎的范围较少，其他领域的专业术语可能翻译不恰当，请各位读者见谅。 **更新状态**： 日期：2022年1月9日 状态：目前全文的内容（除去最后的开源协议）已经通过机翻翻译完毕，后期根据自己需要，会对机翻内容进行删改，形成真正有用的版本。 ## 介绍 ### 什么是Palabos？ Palabos库是一个通用的计算流体动力学 (CFD) 框架，它的内核是格子玻尔兹曼 (LB) 方法。您可以将其当作一个研究工具，也可利用它完成工程项目：其编程界面简单明了，可以相对轻松地进行流体流动模拟；如果您了解格子Boltzmann方法，也可以使用自己的方法对Palabos库进行扩展。 Palabos库是用C++编写的。它几乎没有其他的外部依赖（只有Posix和MPI），因此非常容易部署在各种平台上，包括像BlueGene这样的并行系统。Palabos库还提供了额外的可用于Python和Java的接口，可以快速进行原型设计和开发CFD应用程序。但是，Palabos库目前没有图形用户界面；因此，为了使应用程序运行，我们需要在终端输入一些命令。当使用Python脚本接口时，开发的代码可能非常不正式，但是在我们提供的示例的帮助下，也可以很简单的生成相应的项目。 Palabos库的目标受众包括在计算流体动力学方面具有扎实背景的科学家和工程师，最好具有一些格子Boltzmann方法的一些知识。本软件遵守的开源协议为AGPLv3，其副本印在用户指南的附录中。Palabos库的开源是希望能够促进格子Boltzmann方法领域的研究，帮助研究人员专注于实际的物理问题，而不是停留在繁琐的软件开发中。此外，Palabos库为实现新的格子Boltzmann模型提供了一种简单的方法，以便对新模型进行推广以及交流使用。 Palabos库在高性能计算领域以及复杂流体流动的模拟方面表现出特别出色的性能。在第一种情况下，在拥有数万个内核的机器上使用Palabos库，在测试时我们发现它可以达到接近理想的并行效率。这里一个特别显著的特点是流动模拟的预处理被集成到Palabos的并行管道中。其他软件为一个特定的比较困难的问题去生成网格时可能很困难，但Palabos使用所有可用的计算能力即可在几秒钟内自行构建网格。对于复杂流体流动的模拟，Palabos的多物理场框架可以组合大量的物理模型，同时保持原有的并行效率。使用Palabos执行的复杂任务的示例包括热流模拟、通过多孔介质的多相流模拟以及孔隙尺度精度进行模拟。 Palabos的大量开发时间用于制定LB模拟的通用编程概念，该概念在通用性、易用性和数值效率之间提供了适当的平衡。如果考虑到LB方法的核心是由物理因素而不是数值分析的标准所指导，那么这项任务的困难就可以理解了。因此，可以直接为特定物理模型、矩形数值网格和简单边界条件实现LB代码。然而，在流体工程中，为了解决实际的问题，将现有的简单的代码扩展到能够解决到该困难问题的代码时，我们仍需要仔细考虑以及进行大量的编程工作。在开发Palabos时，经常会发现LB代码实现了一个特定的高级功能（例如并行性），但随后发现它们的形状过于僵化，无法允许其他扩展（例如，多相流的两个网格的耦合）。为了解决这个问题，Palabos软件架构的各种成分（物理模型、基础格点、边界条件、几何域、网格之间的耦合、并行性）在很大程度上发展成了正交概念。例如，可以在不了解Palabos的高级功能的情况下制定经典BGK碰撞模型的变体，例如MRT或熵模型，并自动获取代码的并行版本，或是基于新模型的多相流代码。 Palabos的一个中心概念是所谓的“动力学对象”，它与每个流体单元相关联，并决定局部碰撞步骤的性质。粒子间碰撞的局部性是大多数LB模型的关键组成部分，Palabos通过提高此类模型的数值效率承认了这一事实。特定的数据结构也可用于非严格局部的算法，例如多相模型中的特定边界条件或粒子间力。然而，假设这些成分对项目的整体效率模式只有边际影响。因此，虽然它们以合理的效率实现，但被排除在Palabos的低级数字处理策略之外，以支持更高层次的编程构造。在实践中，我们观察到，通过这种方法，我们得到了程序性能上的最小损失，却降低了最终用户程序的可读性。 模拟的网格结构是多块网格。每个块的行为类似于基本的LB实现，而并行性和稀疏域实现等高级功能则通过块之间的特定接口耦合来实现。 Palabos中大量使用了C++的泛型编程技术。泛型编程的目的是提供可用于多种目的的单一代码。一方面，代码通过面向对象的机制实现了动态泛型。这种机制的一个例子是由格点提供的，它们在程序执行期间改变它们的行为，例如区分体积单元和边界单元，或动态修改流体粘度或体力的值。另一方面，C++的模板用于实现静态泛型。因此，编写一个三维通用代码就足够了，然后在D3Q15、D3Q19和D3Q27晶格上运行。 ### Palabos涵盖的功能 #### 目前已经实现的功能 目前发布的Palabos版本已经完成了以下功能： **物理**：不可压缩NS方程，弱可压缩、等温NS方程，含有体积力的流动，Boussinesq近似的热流模拟，单组分多相流（Shan/Chen模型），多组分多相流（Shan/Chen模型、He/Lee模型），自由表面流动（流体体积方法），湍流中的Smagorinsky模型。 **基本的流体模型**：BGK模型（以及对应的不可压缩变体），MRT模型（以及对应的不可压缩变体），正则化模型，LW-ACM模型（以及对应的不可压缩变体），熵模型。 **直边边界条件**：Zou/He边界条件，Inamuro边界条件，Skordos边界条件，正则化边界条件，平衡态边界条件，反弹边界条件，周期性边界条件。以上所有边界条件都适用于具有内/外角的直壁，并可用于实现速度或压力的Dirichlet或Neumann边界条件。反弹边界条件也可用于曲边边界，只是曲边由阶梯状表示。 **Off-lattice边界条件**：GUO模型，广义off-lattice边界条件。体素化的STL文件的自动、大规模并行的off-lattice壁的实例化。 **粒子**：无源标量或相互作用粒子的大规模并行（在并行机器上数十亿个粒子没有问题）模拟。 **网格**：实现的网格是D2Q9、D3Q13、D3Q15、D3Q19 和 D3Q27。在所有情况下，计算域要么是规则的，要么是一个稀疏域，由多重网格模式近似。 **并行**：所有提到的模型和内容都使用MPI并行化，用于共享内存和分布式内存平台，包括根据MPI的并行I/O接口实现的I/O操作。 **预处理**：模拟域可以手动构建，也可以从相应的STL文件自动构建。 **后处理**：Palabos能够将数据保存在ASCII或二进制文件中，或直接生成GIF图像。此外，数据可以保存为VTK格式，并使用适当的工具进一步进行后处理。为了提高效率，Palabos可以对数据进行本地后处理，生成流线和等值面。 **检查点**：每时每刻都可以保存模拟的状态，并在稍后加载。 #### 正在开发的功能 1. 网格加密。 2. 使用扩展到邻域的晶格模拟热流。 ### 项目 [日内瓦大学](https://www.unige.ch/)负责管理Palabos代码的开发。它仔细审查每一行贡献的代码，判断实现的模型和算法的质量，并促进代码的工业和工程使用。作为计算流体动力学和格子 Boltzmann建模领域的先驱，日内瓦大学的[科学与并行计算小组SPC](http://spc.unige.ch/)提供了Palabos内核的理论背景，并在复杂流体动力学领域不断开发新的模型和方法。 ### 作者 项目监督和核心开发由Jonas Latt负责。贡献者有Orestis Malaspinas， Dimitrios Kontaxakis，Andrea Parmigiani，Daniel Lagrava，Yann Thorimbert，Christos Kotsalos，Francesco Marson，Joel Beny以及Bastien Chopard等。 ### 获取帮助 在开始使用Palabos小节中描述了如何安装Palabos以及代码的简单概述。 除了本文的用户指导手册外，也可以从以下几个地方获取帮助： [教程](https://palabos.unige.ch/get-started/palabos-tutorial/)：为了使用Palabos，强烈建议从教程1开始。该教程提供了高级结构的概述以及有关Palabos内部工作方式的信息。 [论坛](https://palabos-forum.unige.ch/)：论坛可以发布在本用户指南中未回答的问题，也可以与其他Palabos用户、Palabos作者进行交流。 [自动代码文档](https://unigespc.gitlab.io/palabos/docs/html/):该指南是使用[Doxygen](http://www.doxygen.org/)程序从源代码自动生成的。它几乎没有教学价值，但提供了对类层次结构和函数调用语法的详细洞察。 ## 开始使用Palabos ### 支持的编译器 我们没有提供可以成功编译Palabos的编译器和编译器版本号的详细列表。但是，我们已经测试了GCC、Intel编译器和Portland Group编译器，如果您使用最新的版本，它们应该都可以很好地工作。 ### 在Linux或其他Unix系统下安装编译Palabos 最新版本的Palabos可以通过[https://palabos.unige.ch/](https://palabos.unige.ch/)下载。下面我们将描述如何在Linux或Unix系统下如何编译运行Palabos。Palabos库免除了外部依赖：你所需要的只是一个可工作的编译器环境，包括一个现代的C++编译器(例如免费的GCC编译器)，make编译指令以及一个最新版本的Python解释器。首先，下载Palabos最新版本的压缩包，并将其**解压**到您希望放置目录中(我们将以`$(PALABOS)`的方式引用此目录)。 Palabos库使用按需编译过程。代码在用户第一次使用程序时进行编译，之后将自动重用，直到由于修改代码或编译选项而需要进行新的编译。为了展示这一过程，请移动到任意一个示例程序目录下，如`$(PALABOS)/examples/showCases/cylinder2D/`。输入命令`make`，会自动使用GCC编译器编译Palabos库以及对应的项目代码，之后输入`./cavity2D`命令运行代码。该程序模拟了二维圆柱绕流；如果系统中安装了ImageMagick软件，程序会将GIF图像保存在`tmp/`子目录中，表示所选时间步长的速度场。 如果您想要使用不同的编译器进行编译、并行执行，请修改编译选项，只需在本地目录中编辑Makefile文件。以下选项经常被修改: **debug**：打开此标志可以在调试模式下编译。生成的可执行文件会稍微慢一些，但是更容易查找bug。默认打开。 **MPIparallel**：打开并行，最终程序会并行执行。 **serialCXX**：指定串行编译的编译器。 **parallelCXX**：指定并行编译的编译器。 **optimFlags**：编译器优化选项，即`-O0,-O1,-O2,-O3`。 **compileFlags**：指定在特定硬件环境其他编译器选项，例如在Blue-Gene或Mac上编译时，典型的选项分别是-DPLB_BGP或-DPLB_MAC_OS_X。 如果您的系统上安装了MPI（如果您有几个可用的内核或处理器），那么您可以尝试并行编译代码，并通过`mpirun -np 2 ./cavity2d`类型的命令（这行代码表示使用两个线程）执行它。请注意，示例程序的执行时间由输出操作控制。为了观察并行执行中速度的显著提高，首先需要在源代码中将输出注释掉。 ### 在Mac OS X系统上编译Palabos 在Mac OS X系统上，我们可以很容易地通过[xcode](https://developer.apple.com/xcode/)获得GCC编译器，之后的步骤就类似于在Linux系统下的步骤。当然，您也可以仿照在Windows系统下的步骤完成。**提醒**：不管是按照哪种方式，都需要修改Makefile文件，即在Makefile文件中需要使`compileFlags = -DPLB_MAC_OS_X`。 ### 在Windows系统下（或者其他操作系统）使用Code::Blocks集成编译环境编译Palabos C++的集成编译环境（IDE）Code::Blocks是免费的，并且支持多平台使用，具体可见[codeblocks](http://www.codeblocks.org/)。随Palabos一起发布的有一个配置文件，它可以用来通过Code::Blocks编译Palabos。这种方法已经在Windows和Linux下进行了测试，编译过程不需要Python和make。 在Code::Blocks中打开项目文件`$(PALABOS)/codeblocks/Palabos.cbp`，编译和运行这个项目。默认情况下，`examples/showCases/cavity2d.cpp`会被编译。如果您想要编译其他程序，从源文件中移除`cavity2d.cpp`文件，并将您自己的代码文件添加进入即可。生成的GIF图像以及其他输出文件放置在`$(PALABOS)/codeblocks/tmp/`目录下。 Code::Blocks是一个交叉编译器环境，原则上你可以在Code::Blocks中使用任何C++编译器。您可以在Code::Blocks中使用免费的GCC编译器。如果你下载的Code::Blocks版本直接打包了包含GCC编译器的MinGW，您就不需要在Windows系统下单独安装GCC编译器了。 但是请注意，Palabos代码不能用Visual C++编译，我们也不知道为什么，如果您有想法，欢迎向我们提出任何建议。同时，您也可以使用另一种编译器，例如GCC、Intel编译器或Portland Group编译器。 综上所述，Palabos的发行版提供了可以使用免费的集成编译环境Code::Blocks打开的项目文件，您可以使用任何您喜欢的编译器（但是我们只测试了MinGW下的GCC）。如果您的系统上没有安装C++编译器(除了Visual C++)，请您下载将MinGW打包好的Code::Blocks版本，之后您就可以直接编译Palabos了。 如果您希望生成GIF图像的话，请您下载安装[ImageMagick](https://imagemagick.org/index.php)软件。如果您没有安装，Palabos会将图片输出为PPM格式，这个格式在Windows系统下不能被主流图像识别软件识别。在这种情况下，您可能需要安装[Gimp](https://www.gimp.org/)，这是一款图像处理软件，可以查看PPM格式的图片并将其转换为其他格式。 还需要提醒的是，Palabos示例中的路径名是使用Unix约定编写的，在目录名之间使用斜杠，而Windows中使用反斜杠。在Windows系统下，建议这样更改： ``` global::directories().setOutputDir("./tmp/"); ``` 更改为 ``` global::directories().setOutputDir(".\\tmp\\"); ``` （虽然如果您不这样做，Palabos也能正常运行。） ### 在BlueGene/P中编译执行Palabos 在BlueGene/P系统下，编译过程和在Linux或Unix系统下一样。但是要记得修改Makefile文件，设置`compileFlags = -DPLB_BGP`。在BlueGene系统下，GCC的编译速度和XL编译器速度一样快（甚至更快）。 必须指出的是，BlueGene/P系统使用大端法表示数值，而x86架构使用小端法表示数值。因此，您不能将二进制检查点文件从BlueGene导出到x86计算机。BlueGene系统下我们使用另一种方式生成VTK文件，保证能在Intel的PC下查看：只需在您的VTK文件中将字符串`LittleEndian`替换为`BigEndian`即可。您可以通过下面这种方式实现：`sed -i "s/LittleEndian/BigEndian/g" myOutputFile.vti`。对于读取的STL网格，最简单的方法是将STL文件以ASCII格式提供，这种格式与平台无关。如果您的STL文件格式是二进制，您可以通过位于`utility/stl`目录下的`toAsciiSTL`程序将其转化为ASCII格式。 ### 与Palabos绑定的开源库 Palabos利用了其他一些开源库。您不需要手动安装它们，因为它们是Palabos发行版的一部分: [SConstruct](https://www.scons.org/)(仅在Linux系统下进行了测试)：它是一个基于python的库，在Palabos中被用来管理编译过程。 [TinyXML](http://www.grinninglizard.com/tinyxml/)：这个库用于读取XML格式的用户输入。 ### 一些推荐的开源软件 尽管Palabos库是独立的，但是将它与其他产品结合起来是有意义的，特别是对于数据的预处理和后处理。下面列出了一些推荐的开源程序: [Paraview](https://www.paraview.org/)(Linux / Mac OS X / Windows)：Paraview为VTK库提供了一个图形界面，并且是可视化科学数据的优秀工具。Palabos仿真结果可以保存为VTK格式，然后用Paraview进行后处理。在一些Linux发行版的软件存储库中可以找到Paraview。 [Octave](http://www.gnu.org/software/octave/)(Linux / Mac OS X / Windows)：这个命令行解释器与Matlab非常相似。它对于处理Palabos生成的ASCII数据非常有用，可以创建示例图。在大多数Linux发行版的软件存储库中都可以找到Octave。 [ImageMagick](https://imagemagick.org/index.php)(Linux / Mac OS X / Windows)：不需要使用外部库，Palabos可以生成PPM格式的快照图像。如果安装了ImageMagick，这些图片可以自动转换成更常见的GIF格式。此外，还可以使用ImageMagick从一系列GIF图像生成一个动画GIF。 [Meshlab](https://www.meshlab.net/)(Linux / Mac Os X / Windows)：Meshlab可以将许多不同格式的表面网格文件转换为STL格式，这是Palabos设置几何图形所需要的。此外，Meshlab经常能够纠正STL文件中的错误，例如，删除零面积三角形或恢复表面法线。 ### 示例程序 Palabos的`examples`目录分为三个部分：`showCases`，`codeByTopic`以及`tutorial`。`showCases`目录包含了简单流动的全功能实现，`codeByTopic`目录包含了更多的技术代码，以说明在Palabos中特定的编程概念。`tutorial`目录中的内容单独在[palabos_tutorial.pdf](./palabos_tutorial.pdf)文件中描述。 示例程序在附录A中列出和说明。 ## Palabos的Python接口 重要提示:要启动并运行Palabos，不需要编译Palabos的Python接口。Python接口是一个有用的扩展，您会发现它对更快地开发Palabos应用程序很有用。 自从0.7版本以来，Palabos提供了一个Python接口，提供了对库的简单脚本访问。Python接口的外观和感觉与C++接口非常不同，在将来可能会为这个接口提供一个单独的用户指南。同时，本节提供一些编译和使用Python接口的说明。 ### 编译Python接口 为了使用Palabos的Python接口，您需要安装一些库，之后编译Palabos去创建接口。所需的库是所有主要Linux发行版的一部分，可以通过发行版的包轻松安装。当您的包管理器提供了普通版本和`-dev`版本的库时，你必须同时安装它们。需要的库或程序有： * Swig * NumPy * SciPy * Matplotlib * Mpi4py（您可能需要使用`easy_install`安装它。在Ubuntu系统下，通过安装`python-setuptools`可以获得`easy_install`，之后输入`sudo easy_install mpi4py`即可。） 当然了，您还需要安装Python。如果您想要获得更棒的图像，您可能还需要安装`Mayavi2`。 要编译Python接口，您需要一个可工作的MPI和一个可用于MPI的C++包装器(这通常通过包管理器安装`mpich`或`openmpi`来完成)。之后，按顺序完成以下步骤： 1. 找到Palabos的放置位置。设置一个环境变量`PALABOS_ROOT`，存储这个位置。在bash下，需要输入类似`export PALABOS_ROOT=/home/joe/palabos-0.7v0/`的命令。 2. 进入`pythonic/src`子目录下，输入`make`。 3. 如果编译成功，进入`pythonic/examples`目录下，之后输入您想要执行的一个程序即可。 注意，这个编译时间很长，可能需要半个小时。如果出现`file-not-found`的错误，可能是因为没有正确指定库的安装位置导致的。您可以通过找到问题文件所在的目录，并将该目录添加到包含路径来解决这个问题。为此，您需要创建一个Makefile文件，正确编辑`includePaths=...`。 ### 在Ubuntu系统下编译Python接口 以下一步一步的指导是在Ubuntu的9版本和10版本下完成的。下文中出现的具体名字可能需要按照您自己安装的版本进行调整。 #### 安装需要的包 打开包管理器，需要安装下列包：g++，mpi-default-bin，mpi-default-dev，swig，python，python-dev，python-setuptools，python-numpy，python-matplotlib以及mayavi2. 之后打开一个终端，输入`sudo easy_install mpi4py`安装Mpi4Py库。 #### 编译Python接口 打开一个终端，解压下载的压缩包，例如`tar -xvfz palabos-0.7v0.tgz`，进入python接口目录下，例如`cd /home/joe/palabos-0.7v0/pythonic/src`。设置环境变量，即`export PALABOS_ROOT=/home/joe/palabos-0.7v0/`。输入`make`进行编译。 #### 执行示例程序 进入示例目录，例如`cd /home/joe/palabos-0.7v0/pythonic/exmples`，执行一个程序，如`./cavity2d`。您也可以通过在Python提示符中输入命令以交互方式执行程序。为此，显示其中一个示例的内容，例如`cat cavity2d`，通过输入`Python`命令，运行Python解释器，并手动输入程序的行，在运行过程中检查结果。 ## Palabos的Java接口 重要提示:要启动并运行Palabos，不需要编译Palabos的Java接口。Java接口是一个有用的扩展，您会发现它对更快地开发Palabos应用程序很有用。 自从1.1版本以来，Palabos提供了一个Java接口，提供了对库的简单脚本访问。Java接口的外观和感觉与C++接口非常不同，在将来可能会为这个接口提供一个单独的用户指南。同时，本节提供一些编译和使用Java接口的说明。 ### 介绍 Jlabos是Palabos的一个java包装器，它现在是Palabos的一部分，但编译它并不是必须的。 ### 安装 #### 安装要求 * swig1.3 * javac + java，例如`openjdk-6-jdk` * mpicxx，例如`libopenmpi-dev` * mpirun，例如`openmpi-bin` * python 以下内容是可选的： * octave * subversion * convert * inkscape #### 安装 在此之前，您需要确认安装好Palabos。 新建一个环境变量`PALABOS_ROOT`。例如，您可以将这一行代码添加进`~/.bashrc`:`export PALABOS_ROOT=~/location_of_palabos`。 进入源代码目录下，例如`cd $PALABOS_ROOT/jlabos/src`，编辑Makefile文件，设置`nbOfCore = 4`，使用4个CPU核心去编译。这个数值请不要超过你拥有的核心数量，否则将会变得非常慢。同样在Makefile文件中，设置`swigCompiler = swig`，您当然也可以设置为其他。设置`includePaths = ...`，保证这个路径下，含有`jni.h`和`jni_md.h`文件。 编译Jlabos：进入源代码目录下，例如`cd $PALABOS_ROOT/jlabos/src`，在终端内输入`make`命令开始编译。这将耗费很长时间，因为他也会同时开始编译Palabos。 #### 运行示例代码 进入一个您想要运行的示例代码目录下，输入`make`命令，之后即可进行编译和运行。您可以自定义Makefile文件中的一些元素。但是请不要直接启动SConstruct文件。如果需要为共享库指定另一个位置，请相应地编辑`launch`文件。 #### 查看结果 如果您需要，您可以通过Octave或者Matlab绘制结果，只需输入以下命令：`./plot.m`。Jlabos也可以使用嵌入式库`MLplot`直接生成jpg文件。具体请参阅示例代码。`MLplot`需要来自`ImageMagick`的`convert`程序。 #### 自定义包装 Jlabos是Palabos的一个包装，但是它并不完整。如果您需要更多的功能，您需要查看以下两个目录下的内容： 第一个目录是`$PALABOS_ROOT/jlabos/src/jlabos`，在这个目录下，您可以找到封装到自动生成文件中的高级Java文件；第二个目录是`$PALABOS_ROOT/jlabos/src/swig/pre_compiled`，Swig从`.i`文件中生成的文件会显示在这个目录下。**请不要修改java文件！**您只需要修改`.i`文件，之后重新编译运行即可。但是请注意，想要获得一个很好的结果是很困难的。 如果您想要的功能出现在第二个位置，那么您可以在第一个位置修改文件，以很好地包装自动生成的java文件。 ##### 技术细节 包装器的一部分是使用Swig（Simplified Wrapper and Interface Generator）库生成的。这一步不会产生很有用的java文件。所以我们需要利用Java编写第二个包装器，掩盖第一层的结果，保证用户能有较好的使用体验。Jlabos并不支持Palabos的完整功能，但是会不断改进。 `compilepalabos`目录只是强制编译Palabos的一个虚拟目录，您不要过度关注它。 `swig`目录下含有所有swig需要的`.i`文件。 `plbwrapper`包含了一些使用C++编写的Palabos包装器，这些包装器是由swig生成的。 ## 使用Palabos编写代码 ### 快速概述:编程指南 C++是一种非常自由的语言，它对编程风格和编码结构几乎没有什么限制。因此，仅仅了解C++语言并不能开始使用像Palabos这样的库；您还需要熟悉为该库特定算法选择和编程指南。当您熟悉编程规则之后，您就能更确切的知道每行代码是干什么的，并避免像内存泄漏、不正确的并行、预期的面向对象编程的幕后效果等编程错误。本用户指南尝试使您逐步的了解这些概念。每当引入一个新概念时，需要程序员遵守一定的规则来保证正确的结果。 > 你必须遵守规则，否则你的计算结果是不确定的。 本小节将介绍一些您在使用Palabos时的最基本的编程规则，以澄清当您第一次阅读代码时感到困惑的部分。 #### 基本数据类型 在Palabos中，大多数情况下（有例外情况），我们使用`plint`代替`int`类型，用`pluint`代替`unsigned int`类型。这样做的原因是为了兼容64位的机器，并且很容易理解。32位的有符号整数的范围是`-2147483648~2147483647`，64位有符号整数远超过这个范围。因此，很难使用32位机器去模拟超过这个范围的数据集。与32位和64位相关的错误令人很不愉快，您自然也不想遇到这些问题。在Palabos中，当您的机器是64位时，`plint`代表的是64位有符号整数；当您的机器是32位时，`plint`代表的是32位有符号整数。也就是说，当您使用Palabos编写代码时，尽量使用`plint`和`pluint`，除非您有什么必须用`int`和`unsigned int`的理由。 另外，您必须放弃您经常使用的`cout`输出命令。在MPI并行程序中，每个线程通过`cout`写入终端。为了避免这种情况，并在串行程序和并行程序之间观察一致的行为，请使用`pcout`代替`cout`。类似地，用`pcerr`替换`cerr`，用`pclog`替换`clog`。对于文件输出，将`ofstream`替换为`plb_ofstream`。 #### 内存管理 由于C++没有垃圾收集器，使用操作符`new`动态分配的对象必须使用操作符`delete`手动处理。跟踪这些对象并在适当的时候释放它们是这种语言带来的一个令人不快的负担。与其他库不同，Palabos没有包含外部垃圾收集器或类似的构造来解决这个问题。相反，我们希望您会通过遵守一些规则来避免这类问题。 通常，当Palabos中的函数以指向对象的指针作为参数时，这意味着函数(或对象，在类方法的情况下)承担了对该对象的责任。您不用负责删除它。事实上，您并不允许删除它，或是对它做任何事情（因为您并不知道它们在什么地方被删除）。相反的，当函数的返回值是指针时，您就是它的管理者。您可以对他进行任何您想要的操作，同时记得在恰当的时候删除它。 有许多实例都有围绕着谁删除指向对象的指针的问题，在Palabos中通常采用**一个指针一个对象**的方法来避免。当第二个实例需要访问一个已经在使用的对象时，它会创建自己的副本。为了实现这样的方法，所有的Palabos类(或者至少是那些使用继承的类)都有一个`clone`方法: ```C++ A* object1 = new A; A* object2 = object1->clone(); ``` 每当编写继承自Palabos类型的类时，您都需要使用这个方法。幸运的是，`clone`方法的使用总是类似的。 ```C++ A* A::clone() const { return new A(*this); } ``` 不幸的是C++不能自动处理这种事情，但是相信您不久之后就会养成这种习惯。另外，为了保证代码的一致性和可读性，并不建议通过其他方式定义令一个`clone`方法。有时为了`clone`方法有效，有些类可能还需要您提供一个拷贝构造函数。然而，当您在构建这种类型时，我们假设您知道您自己在干什么，并且您知道您需要定义一个拷贝构造函数。 #### 数组 正确选择用于存储值数组的数据类型取决于数据的大小。大的数据集利用Palabos中的`BlockLattice`类型存储，如格子Boltzmann模拟中的分布函数。这些专门化和并行化的数据类型是您开始使用Palabos的首要原因。实例化这样一个类型的语法如下: ```C++ // Instantiate a nx*ny*nz D3Q19 lattice with double-precision // floating point numbers. MultiBlockLattice3D lattice(nx,ny,nz); // Instantiate a nx*ny*nz matrix of double-precision floating // point scalars. MultiScalarField3D field(nx,ny,nz); ``` 小型的、非并行的数组也是需要的，因为我们需要存储一些模拟的参数或者一些结果。C++的标准模板库`vector`可以动态改变数组大小。`vector`类型对我们而言是很有效的，并且使用起来很简单： ```C++ plint numspecies = 5; // Instantiate a vector of double-precision floating point // numbers with 5 elements. vector viscosities(numSpecies); viscosities[3] = 0.63; ``` 强烈推荐您使用`vector`类型去创建动态数组，而不是使用早期的`new []`操作符，以减少错误，提高可读性。 对于非常小的、非并行的、大小固定的数组，Palabos定义了Array类型。它经常用来存储小的物理矢量，如一个单一的速度： ```C++ // Instantiate a fixed-size array of type double and with // 3 elements. The values can be initialized directly in // the constructor for arrays of size 2 and 3. Array velocity(0., 0., 0.5); velocity[0] = 0.1; ``` 同样的，Array类型也要优于继承自C语言的数组（double velocity[3]）。就效率而言，二者的效率是相同的，因为Array类型是使用模板技术，基于早期的数组形成的。然而他有一些优点，例如，当在调试模式下编译Palabos时，它会执行范围检查，这一特性可以大大减少调试时间。 #### 速度和密度 速度和密度是在格子Boltzmann方法中两个最重要的变量。因此，像`rho`或`u`这样的变量很少出现在Palabos代码中，这是十分令人惊讶(但是，没有什么可以阻止您在最终用户代码中阅读它们)。但是，在Palabos中，您会经常见到`rhoBar`和`j`。变量`j`对于大多数模型而言，都是一阶速度矩（即动量），即`j = rho * u`。`rhoBar`的定义是可以自己决定的。默认情况下，它的定义是`rhoBar = rho - 1`，这么定义主要是为了用于提高该方法的数值精度。实际上，由于密度通常接近于1，所以通过表示一个在0而不是1附近的波动的量，可以在浮点变量的表示中获得有效的数字。 我们不进行深入讨论，这里只给出一些交换`rho/u`和`rhoBar/j`的表示方法： ```C++ // Use custom definitions in lattice descriptor to compute rho from rhoBar rho = Descriptor::fullRho(rhoBar); // Use custom definitions in lattice descriptor to compute rhoBar from rho rhoBar = Descriptor::rhoBar(rho); // Momentum is equal to density times velocity (component-wise) j[iD] = rho * u[iD]; // Velocity is equal to inverse-density times momentum (component-wise) u[iD] = 1./rho * j[iD]; // Inverse-density can be computed right away from rhoBar. Depending on // the custom definitions in the descriptor, this is substantially more // efficient than computing 1./rho, because the Taylor expansion of the // inverse-function is truncated at an appropriate position. u[iD] = Descriptor::invRho(rhoBar) * j[iD]; ``` #### 并行 无论使用的是共享式内存平台还是分布式内存平台，Palabos都使用单程序多数据(SPMD)方法实现并行。这意味着当您将程序编写为串行执行或并行执行时，它看起来是一样的（当然，应该加一个定语“至少”，如果您深究的话）。在分布式内存平台上（或者说“集群”），同一个程序的多个实例在不同的线程中启动。在这种情况下，Palabos的分布式对象(MultiBlockLattice、MultiScalarField和MultiTensorField)是并行化的，也就是说，它们在每个线程中分配的内存量需要除以总线程数。通常，其他数据会在每一个线程上拷贝一份。不过，这条规则也有一些例外，我们在本用户指南中突出显示了这些例外。例如，为了提高效率，约简操作的默认行为(例如计算块格的平均密度)是只将结果存储在主线程中。 ### 使用Palabos进行非介入式程序开发 对于在软件工程方面缺乏经验的物理学家和科学家来说，协作编程通常意味着通过电子邮件交换代码片段，并将它们手工集成到现有代码中。在这种方法中，每个程序员都拥有一个完全独立的代码版本，并调整其组件以适应特定的任务。 本节试图说服您不要在Palabos中采用这种方法，而是像使用一个库一样使用这个项目的代码，而不是根据您的本地需求进行调整的代码模板。有几个原因可以解释为什么这样一种非介入性的工作方式在短期和长期内更有效。首先，Palabos的源代码相当大，有5万多行代码，修改整个代码很快就退化成一项非常耗时的任务。第二个，甚至更有说服力的原因是，一旦你修改了Palabos的本地版本，你的副本就会独立于主版本。随着bug修复、效率改进和新特性的出现，要从Palabos的后续发行版中获利变得非常困难，甚至是不可能。此外，您还失去了与同事轻松交换代码片段的可能性，因为不知怎么地，您不得不一直使用Palabos源代码的完整副本。相比之下，基于官方的Palabos版本只交换两个相对简单的代码文件，其中一个为您自己的LB模型实现的一个新的动态类，另一个则是一个示例应用程序。 假设您正在基于OpenGL库生成高级3D图形。你肯定会考虑如何根据OpenGL提供的接口来表述你的问题，而不是修改OpenGL本身的源代码，对吧?以同样的方式考虑Palabos，你就可以节省大量的时间，即使你被迫遵守为Palabos项目所做的一些决定。 #### 扩展代码但不影响原有代码 假如您想写一个类似于BGK模型的LB模型，它们之间只有很小的区别。您需要做的第一件事就是检查实现的流体模型这一小节，确保它还没有在Palabos中实现。如果真的没有，您可能需要查看`src/basicDynamics/isothermalDynamics.h`和对应的`.hh`代码，然后您会发现，`BGKdynamics`继承自`IsoThermalBulkDynamics`，实现了三个方法，`clone()`，`collide()`和`computeEquilibrium()`。根据前面的讨论，您下一步并不是去修改`BGKdynamics`中的`collide()`方法的实现，因为(还有很多其他原因)您可能会失去在应用程序中同时使用原始BGK动力学模型和您的新模型的能力。相反的，创建两个新的文件，`myNewModel.h`和`myNewModel.hh`，保证包含继承自`IsoThermalBulkDynamics`的`MyNewDynamics`类，之后定义和`BGKdynamics`相同的三个方法，以适应您自己的工作要求。这个实现很简洁：它与Palabos版本的更新并不冲突，并且您可以通过交换两个文件与同事共享它。 您的新文件可以位于任何本地目录，甚至独立于Palabos目录。只要确保在Makefile的关键字includePaths下添加这个目录，保证可以在编译过程中找到它即可。 #### Palabos的哪些部分是可扩展的? 简单的说，这里有三种方法去扩展Palabos： 1. 写一个新的动力学类； 2. 写一个新的数据处理器； 3. 写一个新的格子描述方法。 动力学类用于定义一个新的局部碰撞步骤，而一个新的格子描述方法意味着一个新的格子拓扑(离散速度、权重等)，而数据处理器则用于其他的一切。从技术上讲，您还可以通过编写浮点数的自定义类来调整Palabos的行为，但这是一种不太常见的做法。一旦您编写了新的动力学类、数据处理器或格子描述方法，您当然还希望编写一个简单的代码，例如包装器函数，它可以自动为给定的任务交付到正确的数据处理器，或者自动将数据处理器添加到块格中。作为上面所说的代码的一个示例，`OnLatticeBoundaryConditionXD`类自动将动态对象、数据处理器或两者都添加到块格中，以创建边界条件。 Palabos的所有其他部分都不能扩展。例如，写一个类似于`MultiBlockLatticeXD`的新类，或者重新发明`MultiScalarFieldXD`，都不是一个好主意。当然，没有什么可以阻止您这样做，但是重新定义这样的类打破了Palabos中的模块化的概念。简而言之，数据处理器的自定义实现能够从Palabos未来的更新中获益，比如并行性的新方法，而用户发明的`MultiBlockLatticeXD`版本很难编写，而且后续的Palabos版本也不支持。 好消息是Palabos的可扩展部分提供了很大的灵活性。以应用在单组分多相流和非混相多组分流体的二维和三维的Shan/Chen模型的实现为例。这只是通过编写一个新的格子描述方法和一个新的数据处理器来实现的，而不需要修改Palabos的任何其他部分。而且，对应源代码的大小只占整个项目的1%。 不过，有时候在Palabos中所做的选择似乎有点限制性。例如，当我们用`Boussinesq近似`实现热流体模型时，出现了一个问题。在这种情况下，温度场遵循对流扩散方程，该方程采用线性平衡的LB模型求解。从前面的讨论中，很明显，这只需要定义一个实现对流扩散方程的碰撞步骤的新的动力学类。问题是，分布函数的0阶矩在Palabos中总是被称为“rho”(这是由动力学类的界面施加的)，而温度模型的0阶矩是温度。在这一点上，我们必须打破强硬的规则，在代码的几处地方会出现温度称为“ρ”的情况，有必要时请添加一些注释来确保其他程序员不困惑，最终产品包装成一个简单的函数调用温度温度。编程始终是一种折衷的练习，为了避免重写成千上万行代码值得在一些场合调用温度“rho”这个小异端。 ## 基本数据类型 ### BlockXD型数据结构 `BlockXD`型数据结构主要用于存储模拟变量，对于二维问题来说，我们用`Block2D`型数据结构存储模拟变量，对于三维问题来说，我们用`Block3D`型数据结构存储模拟变量。在接下来的讨论中，我们将使用通用的名称`BlockXD`展开讨论(需要注意的是，在Palabos中实际并不存在`BlockXD`这种类型，存在的是`Block2D`和`Block3D`)。从用户的角度看，一个`BlockXD`结构表示一个常规的2维或3维数据数组(或矩阵)。其实，为了实现这种结构，它们有时被真正实现为常规数组，有时被实现为更复杂的构造，这允许通过稀疏内存实现或基于数据并行模型的并行程序执行来节省内存。 `BlockXD`结构可以用于不同的应用领域。一种是专门用于指定存储在块中的数据类型。为了存储格子Boltzmann方法中的粒子分布函数，以及其他可能的变量(如外力)，您将使用特殊化，其中名称`Block`被`BlockLattice`取代。为了存储空间扩展标量变量，要使用的数据类型是`ScalarFieldXD`的一个变体，而向量或张量值字段存储在`TensorFieldXD`或类似的类型中。`AtomicBlockXD`数据结构本质上代表一个常规的数据数组，而`MultiBlockXD`是一个复杂的结构，其中对应于`BlockXD`的空间部分或全部被较小的`AtomicBlockXD`类型的块覆盖。`MultiBlockXD`和`AtomicBlockXD`具有相同的用户界面，我们建议您在终端用户应用程序中使用更通用的`MultiBlockXD`。对于常规问题，它几乎和`AtomicBlockXD`一样高效，但是它可以用来表示不规则域，并且可以自动并行化。 下图展示了BlockXD的各种特殊化之间的C++继承层次结构: ![Alt BlockXDInheritance](./figs/BlockXDInheritance.jpg) 例如，该图显示了在继承层次结构的每个级别上实现的一些函数。例如，所有“block-lattice”类型的块都有一个方法`collideAndStream()`，无论它们是作为`multi-block`还是`atomic-block`实现的。以同样的方式，所有的`multi-block`都有一个`getComponent()`方法，无论它们是`block-lattice`、标量场还是张量场类型。图中所示的情况被称为“多重继承”，因为最终用户类以两种方式继承BlockXD：一种通过原子与多块路径，一种通过`block-lattice`、标量与张量场路径。注意，多重继承通常被认为是不好的实践，因为它会导致代码出错；然而，到目前为止，我们已经有了积极的经验与上面显示的继承图，因为它很容易使用，并以自然的方式表示一个`BlockXD`的两个方面。 ### 格子描述 所有`BlockXD`构造都是利用了模板技术，根据底层数据类型进行构造的。在实践中，该特性用于在单精度和双精度算法之间切换时，只需要更改终端用户程序中的一个单词即可实现: ```C++ // Construct a 100x100 scalar-field with double-precision floating point values. MultiScalarField2D a(100,100); // Construct a 100x100 scalar-field with single-precision floating point values. MultiScalarField2D b(100,100); ``` `Block-lattices`还有一个模板参数，即格子描述符，它指定了格子的一些拓扑属性(分布函数的数量、离散速度、方向的权重和其他格子常量)。因此，在应用程序中尝试不同的格子是很容易的: ```C++ // Construct a 100x100 block-lattice using the D3Q19 structure. MultiBlockLattice2D lattice1(100,100); // Construct a 100x100 block-lattice using the D3Q27 structure. MultiBlockLattice2D lattice2(100,100); ``` 编写一个新的格子描述符也很容易(该方面的内容目前还没有文档，但您可以查看`src/latticeBoltzmann/nearestNeighborLattices2D.h`目录中的文件，以了解它是如何工作的)。这是非常有用的，因为这意味着当您切换到一种新类型的格子时，您不需要为`BlockLatticeXD`的实现重写冗长的代码部分。这个参数是使用Palabos进行非介入式程序开发一节中描述的一般概念的一部分。 ### 动力学类 在使用格子Boltzmann方法模拟的一个时间步中，所有的格子单元都需要执行一个局部碰撞步骤，之后是一个迁移步骤。迁移步骤是固定的，只能通过定义格子描述符中的离散速度来影响。另一方面，碰撞步则可以完全定制，不同的格子单元可以不同。通过这种方式，在网格上模拟的物理性质可以局部调整，特定的区域如边界可以得到单独的处理。 ![Alt](./figs/InheritanceDiagramForDynamicsClasses.jpg) 要学习如何定义一个新的动态类，最简单的方法是查看在Palabos中定义的一个类。例如，BGK动态类是在文件`src/basicDynamics/isoThermalDynamics.hh`中定义的。复合动力学的一个很好的例子是`Smagorinsky`动力学，它定义在`src/complexDynamics/smagorinskyDynamics.hh`中。 碰撞步的施加，动态对象会使用一个局部格子的引用；因此，碰撞步就可以保证是局部的了。非局部成分的模拟是通过数据处理器实现的，将在下一节描述。 与`block-lattice`类似，一个动态对象依赖于两个模板参数，一个是使用哪种浮点数类型（`float`或`double`），另一个是格子描述符。通过使用格子描述符提供的信息，碰撞步相应的代码应该是通用的，与格子无关。以通用方式编写算法显然存在效率权衡，因为可以为编译器无法编译的特定格子制定优化。这个问题可以通过使用给定的格子模板特化来规避。再以`BGKdynamics`类的实现为例。碰撞步的实现引用了一个通用对象`dynamicsTemplates`，它定义在文件`src/latticeBoltzmann/dynamicsTemplates.h`中。在文件 `dynamicsTemplates2D.h`和 `dynamicsTemplates3D.h`中可以找到此类针对各种2维和3维晶格的有效专业化。 ### 数据处理器 数据处理器定义了要在整个域或块的一部分上执行的操作。在块格上，它们用于实现所有不能用动力学对象表述的操作。这些操作包括最显着的非局部操作，例如评估有限差分模板以实现边界条件。在标量场和张量场上，数据处理器提供了在空间扩展域上执行操作的唯一（足够高效且可并行化）的方法。 此外，数据处理器能够在相同类型或不同类型的多个块之间交换信息（例如：块格和标量场之间的耦合）。例如，这用于实现物理耦合（多组分流体、具有 Boussinesq近似的热流体），用于设置初始条件（从矢量场的值初始化速度），或执行基于经典数组的操作（两个标量字段的元素相加）。 最后，数据处理器用于对整个块或子域执行归约操作。例如计算模拟中的平均动能，计算作用在障碍物上的阻力。 对于数据处理器的定义和应用，采用了两种不同的观点。在应用程序级别，用户指定给定块的一个区域（可以是矩形或不规则的），在该区域上执行数据处理器。如果块内部具有多块结构，则数据处理器被细分为多个更具体的数据处理器，多块内与指定区域相交的每个原子块一个。因此，在执行级别，数据处理器总是作用于原子块，作用于先前通过将原始区域与原子块的域相交而确定的区域。 应该提到的是，虽然原始数据处理器使用起来有些笨拙，但您可能永远不会与它们接触。相反，Palabos通过所谓的数据处理功能提供了一个简化的界面，它隐藏了技术细节，让您专注于基本部分。用户指南的其余部分专门针对这些功能，简称为数据处理器。 最后，定义新的数据处理器非常容易。您需要做的就是编写一个函数，该函数接收原子块和子域的坐标作为参数，并在该子域上执行算法。所有复杂的操作，比如存在多块时的操作细分，或代码的并行化，都是自动的。 在示例文件`examples/codeByTopics/couplings`中可以找到数据处理器的示例。它展示了如何通过为块点阵与二维张量场耦合编写数据处理器来使用矢量场中的速度初始化块点阵。 作用于块格的数据处理器的更多定义可以在文件`src/simulationSetup/latticeInitializerXD.h` 和对应的`.hh`中找到，作用于标量或张量场的数据处理器在文件`src/simulationSetup/dataFieldInitializerXD.h`及相应的`.hh`文件中定义。 文件`src/core/dataAnalysisXD.h`和对应的`.hh`文件中提供了用于评估归约操作的示例。 所有这些数据处理器都包含在方便的函数中，这些函数在附录中进行了总结，请参考部分函数/类部分。 ## 实现的流体模型 ### 等温NS方程 有一些解决可压缩流体动力学的格子Boltzmann模型，但由于各种原因，这些模型中关于能量方程的建模是错误的。这些模型有时会用于温度变化可忽略不计的小马赫数可压缩流动，但大多数情况下，它们仅用于求解不可压缩的Navier-Stokes（NS）方程。在这种情况下，马赫数保持较低以抑制不需要的可压缩性效应，并且该模型用作所谓的准可压缩求解器。 除非另有说明，否则所有呈现的模型都可以与D2Q9、D3Q15、D3Q19和D3Q27离散速度模型一起使用，我们使用通用的描述符`DdQqDescriptor`描述这些离散速度模型。 如果你使用了含有外力项的模型，请用`ForcedDdQqDescriptor`替换`DdQqDescriptor`。 应该提到的是，在Palabos开发人员中，熵模型存在不确定性（它可能仅在D2Q9和D3Q27上完全有效）。 #### 单松弛BGK模型 这个模型是最基础的LB模型，它被广泛使用，实施起来很简单简单并且用途很广泛。 | Dynamics for plain BGK: | BGKdynamics | |:-----------------------------:|:--------------------------------:| | Dynamics with external force: | GuoExternalForceBGKdynamics | | | ShanChenExternalForceBGKdynamics | | | HeExternalForceBGKdynamics | | Example: | codesByTopic/navierStokesModels | #### 不可压BGK模型 在标准BGK模型中，平衡态项需要乘以流体密度。在所谓的不可压缩模型中，密度仅出现在常数项中（该项不含流体速度），速度正比于动量大小（注意不需要除以密度）。这个模型相比BGK模型来说，计算起来更加简单方便，因为在计算速度时不需要除以密度。此外，与不可压缩NS方程的精确解相比，模拟的可压缩性误差对于定常流减少了，因为它随着马赫数的立方一样减小，而不是标准BGK模型中的马赫数平方。但是，在非定常流动中，误差仍然与马赫数平方成正比。 | Dynamics: | IncBGKdynamics | |:---------:|:-------------------------------:| | Example: | codesByTopic/navierStokesModels | #### 正则化BGK模型 在碰撞步之前，粒子分布函数被投影到具有二阶多项式的Hermite基上。虽然这保持了NS方程的动力学完整，但它通过对粒子分布函数施加额外的稳定性来提高数值稳定性。 `CompositeDynamics`类型的正则化碰撞存在一个通用版本，可用于正则化任何LB模型。此外，出于效率原因，为正则化BGK模型实现了一个特定的类。 | Generic Dynamics: | RLBdynamics | |:-----------------:|:-------------------------------:| | Dynamics for BGK: | RegularizedBGKdynamics | | Example: | codesByTopic/navierStokesModels | #### 线性稳定性优化的多松弛模型（MRT） 在多松弛（MRT）模型中，线性碰撞算子被重铸到速度矩空间中，并单独调整它们相应的松弛参数，以修改模型的物理特性，或提高数值稳定性。在这里实现的模型中，特定参数用于非物理重影模式以提高稳定性。 | Dynamics: | MRTdynamics | |:------------:|:-------------------------------------:| | Descriptors: | MRTD2Q9Descriptor, MRTD3Q19Descriptor | | Example: | codesByTopic/navierStokesModels | #### 熵模型 在没有边界的情况下，这个模型是无条件稳定的，因为它保证了熵函数的存在，就像自然界一样，它是非递减的。 | Dynamics for plain model: | EntropicDynamics | |:-----------------------------:|:-------------------------------:| | Dynamics with external force: | ForcedEntropicDynamics | | Example: | codesByTopic/navierStokesModels | ### 使用Boussinesq近似的热流 在`Boussinesq`近似中，流体被认为是不可压缩的，温度的影响只能通过代表浮力效应的体积力项可见。温度场服从对流扩散方程。 温度场和流体求解器之间的耦合是通过数据处理器实现的。对于流体，您可以使用任何具有外力项的等温NS模型。 #### 温度的对流扩散模型 | Dynamics: | AdvectionDiffusionBGKdynamics | |:-----------------------------:|:-------------------------------:| | Descriptors: | AdvectionDiffusionD2Q5Descriptor, AdvectionDiffusionD3Q7Descriptor | #### 流体与温度的耦合 | Data Processor: | BoussinesqThermalProcessorXD | |:-----------------------------:|:-------------------------------:| | Examples: | showCases/boussinesqThermal2D 以及 showCases/boussinesqThermal3D | ### 多组分多相流 广泛用于不混溶多组分流体和单组分多相流体的Shan/Chen模型已经实现了在Palabos中的2维和3维问题的应用。请注意，有关这些模型的理论和实践方面的大量信息可以在[Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers](https://link.springer.com/book/10.1007/978-3-540-27982-2)书中找到。 在这些模型中，都需要计算一个并非完全局部的粒子间作用力。因此，该作用力项由数据处理器计算。该数据处理器需要知道给定格点上的密度和速度（在多分量模型中，必须知道所有分量的这些宏观变量）。为了避免对这些变量进行两次计算，因为它们也需要在后续的正常碰撞步骤中知道，数据处理器将它们写入单元的外部标量，然后由动力学对象重新使用以进行碰撞。因此，有必要使用为这些外部标量提供空间的格子描述符，例如`ShanChenD2Q9Descriptor`或`ForcedShanChenD2Q9Descriptor`，如果除了粒子间势能之外还有外部体力。在3维问题中，相同的描述符可用于D3Q19点阵（并且其他描述符很容易制定）。为确保在碰撞步期间从外部标量读取宏观变量，请为流体模型使用`ExternalMomentBGKdynamics`或`ExternalMomentRegularizedBGKdynamics` 类。 请注意，在两种模型中，如果在格子描述符中定义了外力，则会自动考虑外力项的影响。然后，Shan/Chen模型对应的数据处理器通过对流体速度进行动量校正来添加力项。因此，您不需要使用像`GuoExternalForceBGKdynamics`这样的带有力项的特定碰撞模型。简单的`BGK`模型或正则化的`BGK`模型就已经足够了。 #### 多组分Shan/Chen模型 在该模型中，每个组件都在单独的块格子上进行仿真，并且这些组件通过Shan/Chen数据处理器进行耦合。您可以根据需要连接任意数量的组件。 | Data Processor: | ShanChenMultiComponentProcessorXD | |:---------------:|:---------------------------------------------------------:| | Examples: | showCases/multiComponent2d 以及 showCases/multiComponent3d | #### 单组分Shan/Chen多相模型 在这个模型中，Shan/Chen数据处理器采用了一个额外的参数来定义粒子间势的形状。定义了以下电位（有关详细信息，请参阅文件`src/multiPhysics/interparticlePotential.h`）： | PsiIsRho() | $\Psi=\rho$ | |:------------------------:|:-:| | PsiShanChen93(rho0) | $\Psi=\rho_{0}\left( 1-e^{-\rho/\rho_{0}} \right)$ | | PsiShanChen94(Psi0,rho0) | $\Psi=\Psi_{0}e^{-\rho_{0}/\rho}$ | | PsiQian95(rho0,g) | $\Psi=g\rho_{0}^{2} \rho^{2} / (2(\rho_{0} + \rho)^2)$ | 请注意，在此模型中，流体密度必须调整以适应相互作用力的幅度，以便进入发生相分离的临界状态。 | Data Processor: | ShanChenSingleComponentProcessorXD | |:---------------:|:----------------------------------:| | Examples: | codesByTopic/shanChenMultiPhase | ### 自由表面流动 在自由表面流动中，两相之一的粘度实际上是无穷大。因此该阶段被忽略，而表面跟踪方法被应用于剩余阶段。Palabos使用类似流体体积的方法来追踪自由表面。 例如，在`examples/showCases/breakdownDam3d`中提供了一个破坝应用程序。 ### 大涡模拟 #### 静态Smagorinsky模型 在Smagorinsky模型中，假设亚网格尺度具有粘度校正的效果，该效果与过滤尺度级别的应变率张量的范数成正比，即$\nu = \nu_0 + \nu_T$。涡粘性系数计算公式为 $$ \nu_{T} = C^2|S| $$ 其中，$C$是Smagorinsky常数，应变率的张量范数定义为$|S| = \sqrt{S:S}$。Smagorinsky常数的值取决于我们所模拟的物理问题，远离边界的区域通常的取值范围在0.1到0.2之间。这个模型是一个静态模型，因为Smagorinsky常数是强加的并且不会随时间改变。 在Palabos的实现中，应变率$S$是从应变张量$\Pi$计算得来的。值得注意的是，应变率$S$和应变张量$\Pi$之间的关系和松弛时间$\tau$有关，从而与运动粘性系数$\nu$有关。涡粘性系数$\nu_{T}$的公式是隐式的，但是由于它是二阶的，是可以明确求解的。 鉴于可以从单元上的局部变量计算应力张量$\Pi$，Smagorinsky模型完全是局部的，并通过动力学类在Palabos中实现。通常，此模型有多种实现方式。 通用的基于复合动力学对象，在运行时向任何现有动力学类中添加了Smagorinsky粘度校正。特定的那些是为BGK和正则化BGK模型明确编写的，因此效率更高。在每种情况下，在执行通常的碰撞步之前修改粘度值。因此，如果在模拟过程中您访问单元中的松弛参数$\omega$，例如通过调用`lattice.get(x,y).getOmega()`的函数，您将获得与有效粘度$\nu$相关的松弛参数，而不是 “分子粘度”$\nu_{0}$。 Smagorinsky常量的值必须作为构造函数参数提供给这些动态类。需要注意的是，每个单元的动力学对象可以具有不同的Smagorinsky常数值：这个“常数”可以是空间相关的。例如，这对于模拟边界层很有用，其中Smagorinsky常数降至零。创建具有空间相关Smagorinsky常数的模拟的最方便的方法是从数据处理器中将 Smagorinsky 动力学分配给晶格的每个单元格，该处理器知道作为空间函数的$C$的所需值。 | Generic Dynamics: | SmagorinskyDynamics | |:-----------------------------:|:------------------------------:| | Dynamics for BGK: | SmagorinskyBGKdynamics | | Dynamics for Regularized BGK: | SmagorinskyRegularizedDynamics | | Example: | codesByTopic/smagorinskyModel | ### 非牛顿流体 #### Carreau模型 在Carreau模型中，粘性的值是可调整的，与Smagorinsky模型类似，粘性的大小依赖于当地的应变率张量。本构方程由下式给出 $$ \nu=\nu_{0} \times \left(1+(\lambda \times |S|)^{2}\right)^{(n-1) / 2} $$ 其中，$\lambda$和$n$是给定的实参。在Smagorinsky模型中，应变率张量$S$是从应变张量$\Pi$中获得的，因此最终的涡粘性系数的计算是隐式的，因为应变率张量$S$和应变张量$\Pi$的关系取决于松弛时间$\tau$的大小。然而，在目前的情况下，并没有明确的解，方程在每个时间步需要通过点松弛迭代求解（在这种情况下收敛非常快，因此比基于梯度的求解方法更有效）。 在具体实现中，$\nu_{0}$，$\lambda$和$n$等参数并不是空间相关的，它们是通过对全局单例`CarreauParameters`的函数调用来设置的。例如 ```C++ global::CarreauParameters().setNu0(nu0); global::CarreauParameters().setLambda(1.); global::CarreauParameters().setExponent(0.3); ``` | Dynamics: | CarreauDynamics | |:---------:|:-------------------------:| | Example: | codesByTopic/nonNewtonian | ## 设定一个具体的问题 ### 分配动力学对象 当构造一个新的块后，您必须通过动力学对象确定在这一区域的每一个单元要进行的碰撞类型是什么样的，主要是确定这一区域每一个单元的松弛频率是多少。为了避免出现由于您没有动力学对象导致的Bug，每一块在构造时会默认附带一个动力学对象，被称作“默认动力学类”： ```C++ // Construct a new block-lattice with a background-dynamics of type BGK. MultiBlockLattice3D lattice( nx, ny, nz, new BGKdynamics(omega) ); ``` 在这之后，每一个单元的动力学对象都可以重新定义，以保证符合您自己的要求。默认动力学类与一般的动力学类有很大的不同。为了节约内存，块在构造时，只会创建一个动力学对象，所有的单元都指向这一个动力学对象。如果你在一个单元上修改了动力学对象，整个块都会被修改。以松弛频率$\omega$为例，它被存储在动力学对象中，并不是在每一个单元中。一个函数调用如下 ```C++ lattice.get(0,0,0).getDynamics().setOmega(newOmega); ``` 该代码会影响整个块中所有单元的松弛频率。注意，这是针对默认动力学对象的特殊现象。除非您创建一个新的块，否则这种现象是不可能重新出现的，即不可能出现多个单元引用了同一个对象。使用新动力学对象覆盖给定单元格上的默认动力学对象的所有其他函数的定义方式是为每个单元创建动力学对象的独立副本。 如果默认动力学的引用语义与您的需要相反，您可以在创建块后立即重新定义所有动力学对象，以保证所有单元都有一个独立的副本： ```C++ // Override background-dynamics to guarantee and independent per-cell // copy of the dynamics object. defineDynamics( lattice, lattice.getBoundingBox(), new BGKdynamics ); ``` 函数`defineDynamics`存在不同的变体，您可以使用它们来调整一组单元的动力学对象（它们的语法可以在附录块的操作中找到）： * **One-cell版本**：仅将新的动力学对象与一个单元关联。示例代码：`tutorial/tutorial2/tutorial2_3.cpp` * **BoxXD版本**：将新的动力学对象与一个Box区域内的所有单元关联。示例代码：`showCases/multiComponent2d/rayleighTaylor2D.cpp`，或与其对应的三维代码 * **DotListXD版本**：将新的动力学对象与多个单元关联，这些单元以点列表结构单独列出。示例代码：`codesByTopic/dotList/cylinder2d.cpp` * **Domain-functional版本**：提供一个分析函数，该函数指示获得新的动力学对象的单元的坐标。示例代码：`showCases/cylinder2d/cylinder2d.cpp` * **Bool-mask版本**：通过布尔掩码指定获得新的动力学对象的单元的位置，通过标量场表示。示例代码：`codesByTopic/io/loadGeometry.cpp` ### 初始化速度和密度 通过将所有单元初始化为给定密度和速度的平衡态分布函数，将速度和密度的初始值分配给每一个单元是很常见的。虽然这种方法对于严重取决于初始状态的依赖于时间的问题是不够的，但通常足以以合理的初始值开始模拟。提供了`initializeAtEquilibrium`函数（请参阅附录块的操作）以在均衡分布的矩形域内初始化单元。它有两种风格：一种在域内具有恒定的密度和速度（参见示例`examples/showCases/cavity2d`或对应的`3d`代码。），另一种具有这些宏观量的空间相关值，通过用户定义的函数指定（参见示例`examples/showCases/poiseuille`）。 要以不同的方式初始化单元格，最简单的方法是在单元功能的帮助下将自定义运算符应用于单元格，这些功能在一些对于本地操作很便利的包装器中介绍。 ## 设定边界条件 ### 概览 Palabos库目前针对直边、网格对齐边界问题已经实现了很多种边界条件。对于一般的边界，可用选项包括通过反弹节点的楼梯墙，以及高阶精确曲线边界。 重要提醒：曲边边界从1.0版开始即可用，并且是功能齐全的并行堆栈的一部分，包括并行体素化和I/O扩展。但是曲边边界暂时没有说明文档。您可以在`showCases/aneurysm`中找到功能齐全的示例。 ### 网格对齐边界 #### `OnLatticeBoundaryConditionXD`类 一些边界条件的施加是局部的（因此一般按照动力学类的方式施加），一些则是非局部的（因此一般按照数据处理器的方式施加），还有一些是二者混合的。针对这些情况，为了提供一个通用的接口，Palabos提供了`OnLatticeBoundaryConditionXD`接口，它负责实例化动力学对象，添加数据处理器，或两者兼而有之，具体取决于所选的边界条件。下面是使用Zou/He边界的一个示例： ```C++ OnLatticeBoundaryCondition3D* boundaryCondition = createZouHeBoundaryCondition3D(); ``` 目前我们提供的有： | Regularized BC | createLocalBoundaryConditionXD | | :--------------: | :-------------: | | Skordos BC | createInterpBoundaryConditionXD | | Zou/He BC | createZouHeBoundaryConditionXD | | Inamuro BC | createInamuroBoundaryConditionXD | 对于这些边界条件中的每一个，都可以实现速度Dirichlet边界（所有速度分量都有一个强加值）和压力边界（强加压力，切向速度分量为零）。另外，许多不同类型的Neumann边界也被提出。在这些情况下，通过从相邻单元复制该变量的值，给定变量的零梯度以一阶精度强加。 #### 规则区域的边界 这些边界条件的使用有些棘手，因为Palabos需要知道给定的壁面是直壁面还是拐角（2D），还是直壁面、边或拐角（3D），它需要知道壁面方向。为了简化这种情况，所有的`OnLatticeBoundaryConditionXD`对象都有两个方法去设置位于规则域上的边界，分别是`setVelocityConditionOnBlockBoundaries`和`setPressureConditionOnBlockBoundaries`。如果给定二维框`boundaryBox`边界上的所有节点都施加Dirichlet条件，请使用以下命令： ```C++ Box2D boundaryBox(x0,x1, y0,y1); boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, boundaryBox, locationOfBoundaryNodes ); ``` 第一个参数表示边界的形状，第二个参数指定要在哪个节点上施加速度条件。如果所有边界节点都位于计算域的外部边界框上，则上述命令简单地变为： ```C++ boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes); ``` 最后，如果简单地计算域的所有节点都应实现速度条件，则简化为： ```C++ boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice); ``` 现在，假设 $nx \times ny$ 计算域中的上壁面和下壁面（包括角）实现自由滑移条件（切向速度分量为零梯度），左壁面是Dirichlet条件，右壁面是流出条件（所有速度分量的零梯度）。这是通过在函数调用中再添加一个参数来实现的： ```C++ Box2D inlet(0, 0, 2, ny-2); Box2D outlet(nx-1, nx-1, 2, ny-2); Box2D bottomWall(0, nx-1, 0, 0); Box2D topWall(0, nx-1, ny-1, ny-1); boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, inlet ); boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, outlet, boundary::outflow ); boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, bottomWall, boundary::freeslip ); boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, topWall, boundary::freeslip ); ``` 请记住，如果边界不在格子的外边界上，则必须使用额外的`Box3D`参数来说明边界所在的位置，此外还要说明要添加哪个边界。这是必要的，因为Palabos需要知道边界节点的方向和性质（壁面、拐角等）。 #### 零梯度边界条件 `boundary::outflow`是一个可选参数，它是`boundary::BcType`类型，对于速度边界有以下选择： | boundary::dirichlet (default value) | Dirichlet condition: imposed velocity value. | | :--------------: | :-------------: | | boundary::outflow or boundary::neumann | Zero-gradient for all velocity components. | | boundary::freeslip | Zero-gradient for tangential velocity components, zero value for the others. | | boundary::normalOutflow | Zero-gradient for normal velocity components, zero value for the others. | 对于压力边界，您有以下选择： | boundary::dirichlet (default value) | Dirichlet condition: imposed pressure value, zero value for the tangential velocity components. | | :--------------: | :-------------: | | boundary::neumann | Zero-gradient for the pressure, zero value for the tangential velocity components. | --- **边界条件不能被覆盖** 不可能覆盖边界的类型。一旦边界节点是Dirichlet速度节点，它就一直是Dirichlet速度节点。将其重新定义为压力边界节点的结果是未定义的。因此，必须仔细地、分段地定义边界条件，如上例所示。另一方面，施加在边界上的值（即Dirichlet速度边界上的速度值）可以根据需要随时更改。 --- #### 速度和压强的Dirichlet边界条件 `setBoundaryVelocity`函数会强加一个恒定的速度值。下面的命令将二维域的所有节点上的速度值设置为零，这些节点之前已定义为Dirichlet速度节点： ```C++ setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array(0.,0.) ); ``` 在所有其他节点上，此命令无效。要设置恒定压力值（或等效的密度值），请使用命令`setBoundaryDensity`： ```C++ setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. ); ``` 可以通过用函数或函数对象（覆盖函数调用运算符的类的实例）替换速度和压力参数来定义非固定的、空间相关的速度和压力分布，从而产生速度值和密度作为空间位置的函数。文件`examples/showCases/poiseuille/poiseuille.cpp`中提供了一个示例。 #### 内部和外部边界 有时，边界并不位于规则块状域的表面上。在这种情况下，像`setVelocityConditionOnBlockBoundaries`这样的函数就没有用了，必须手动构建边界。例如，以下是如何沿顶壁手动构建自由滑移条件，包括角节点： ```C++ // The automatic creation of a top-wall ... Box2D topWall(0, nx-1, ny-1, ny-1); boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, topWall, boundary::freeslip ); // ... can be replaced by a manual construction as follows: Box2D topLid(1, nx-2, ny-1, ny-1); boundaryCondition->addVelocityBoundary1P( topLid, lattice, boundary::freeslip ); boundaryCondition->addExternalVelocityCornerNP( 0, ny-1, lattice, boundary::freeslip ); boundaryCondition->addExternalVelocityCornerPP( nx-1, ny-1, lattice, boundary::freeslip ); ``` 外角和内角之间有区别，这取决于角是凸面还是凹面。 在以下示例几何体中，您会发现五个外角和一个内角： ![alt](./figs/ConvexorConcave.jpg) 边界条件对象方法末尾的`1P`、`NP`和`PP`等扩展用于指示壁面法线的方向，指向流体域之外，如上图所示。在平直壁面上，代码`1P`表示：“壁面法线指向正y方向”。同样，入口将用代码`0N`标记为“负x方向”。在拐角处，代码`NP`表示“负x方向和正y方向”。值得一提的是，该壁面法线纯粹是几何形状决定的，并不取决于给定的壁面是具有入口还是出口的功能。在这两种情况下，壁面法线都远离流体，指向非流体区域。 在3D中，您将对平面使用以下函数： ```C++ addVelocityBoundaryDO ``` 其中方向D可以是x的0，y的1或z的2，方向O的值P表示正，N表示负。边缘可以是内部的，也可以是外部的。例如： ```C++ addInternalVelocityEdge0NP addExternalVelocityEdge1PN ``` 其中代码的第一位数字表示边平行的轴，后面的两个数字表示与边垂直的平面内边的方向，与二维中的角节点相同。轴会循环计数：0NP表示x平面、负y值和正z值，而1PN表示y平面、正z值和负x值。最后，以与2D相同的方式构建3D角，通过如下所示的函数调用： ```C++ addInternalVelocityCornerPNP addExternalVelocityCornerNNN ``` ### 周期性边界条件 在Palabos中，周期性的概念与其他类型的边界条件不同。周期性边界条件仅能在单块或多块网格的相反外边界上施加。即使多块网格使用了稀疏内存的技术实现，该方法同样可以正常实现。在这种情况下，如果相对壁上的相应节点是流体节点，则与多块的外边界接触的所有流体节点都可以是周期性的。 在单块网格的情况下（请记住，这种情况是无趣的，因为无论如何您都应该始终使用多块网格），周期性只是流操作符的一个属性。它对标量场和张量场没有影响。然而，在多块网格的情况下，周期性也会影响标量场和张量场。在这种情况下，周期性意味着如果您定义访问最近邻居节点的非本地数据处理器，则这些邻居节点的值由沿多块外边界的周期性决定。 默认情况下，所有块都是非周期性的：如前所述，外边界节点的行为默认为反弹算法版本之一。要沿x方向打开周期性，请编写如下命令： ```C++ lattice.periodicity().toggle(0, true); // Periodic block-lattice. scalarField.periodicity().toggle(0, true); // Periodic scalar-field. ``` 您也可以定义所有的边界均为周期性边界条件， ```C++ lattice.periodicity().toggleAll(true); ``` 请注意，您最好在定义完块之后，立刻定义周期性边界条件或者非周期性边界条件，因为Palabos在施加数据处理器的时候需要知道是周期性边界还是非周期性边界。 最后，您应该知道Neumann类型的边界（流出、自由滑移、绝热壁面等）永远不应该定义为周期性的。因为这是没有意义的，它会导致未定义的行为（也就是程序崩溃）。 ### 反弹边界条件 在Palabos中，所有的外部边界，只要您未指定是周期性边界条件，均默认是反弹边界条件。该边界条件施加如下：在碰撞步之前，某一方向未知的分布函数在上一次迭代结束后，均被赋值为同一节点，对应相反方向的碰撞后的分布函数。我们称这种为`Bounce-Back 1`。 反弹条件可以在计算域中的任意位置使用，通过将所选择的单元与一个`BounceBack`类型的动力学对象关联起来即可。在这些新关联起来的节点上，碰撞步被反弹过程替代，即每一个方向的分布函数被其对应相反方向的分布函数替代。这些节点就不能再被视为是流体节点了。计算这些节点上的分布函数的速度矩会产生各种可能，因为有些分布函数（有些分布函数并不是来自流体节点）是任意大小的。所以，您更不能计算这些节点的宏观量，如速度和密度。相反的，您应该将这些节点视为纯粹为了实现相关算法所需要的虚拟节点，这些节点的目的就是为了将所有离开域的种群重新注入流体中。现在我们考虑十分靠近这些反弹节点的流体节点。您会发现这些十分靠近这些反弹节点的流体节点的行为与`Bounce-Back 1`的行为很像，仅有小小的不同：$t$时刻未知的分布函数被赋值为$t-2$时刻碰撞后的分布函数（请注意，`Bounce-Back 1`是$t$时刻未知的分布函数被赋值为$t-1$时刻碰撞后的分布函数）。这时的反弹类型我们称之为`Bounce-Back 2`。 这两种反弹边界条件的最终效果是在两个格点正中间实现无滑移边界条件。在`Bounce-Back 1`中，这是反弹流体单元之外的单元宽度的一半，而在`Bounce-Back 2`的情况下，它位于最后一个流体单元和非流体反弹节点之间的中间。显然，`Bounce-Back 2`通过多跨越一个时间步而降低了`Bounce-Back 1`的数值精度，因此导致该数值方法的时间精度较低。此外，两个版本的反弹都通过楼梯近似表示壁面位置，这是空间中的低阶准确（“Order-h”）。另一方面，`Bounce-Back 2`提供了一个巨大的优势：它的实现与墙壁方向无关。 如果你决定一个节点是一个反弹节点，你只需这么说；无需知道它是在角落、左壁面还是右壁面。构建几何体就像用反弹节点填充区域一样简单。在许多情况下，这种易用性在很大程度上弥补了壁面表示的低阶精度。设置新模拟时，最好先尝试`Bounce-Back 2`，因为所获得结果的质量通常令人惊讶。在高度复杂的几何形状（例如多孔介质）中，反弹边界通常甚至被认为优于其他方法。 #### 分析描述的反弹域 使用Palabos的第一准则是：在您编写的代码中不要出现空间遍历。当然，在您将反弹动力学分配给晶格单元时，您也应该这样。为了保证合理的性能，这是在数据处理器内部执行的任务，或者甚至更好，通过使用包装函数`defineDynamics`。这个过程在示例程序`examples/codesByTopic/bounce Back/instantiateCylinder.cpp`中有说明。首先，编写了一个类，它将反弹节点的位置定义为单元坐标iX和iY的函数。在本例中，节点位于二维圆形域内： ```C++ template class CylinderShapeDomain2D : public DomainFunctional2D { public: CylinderShapeDomain2D(plint cx_, plint cy_, plint radius) : cx(cx_), cy(cy_), radiusSqr(util::sqr(radius)) { } // The function-call operator is overridden to specify the location // of bounce-back nodes. virtual bool operator() (plint iX, plint iY) const { return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiusSqr; } virtual CylinderShapeDomain2D* clone() const { return new CylinderShapeDomain2D(*this); } private: plint cx, cy; plint radiusSqr; }; ``` 然后将此类的实例作为参数提供给`defineDynamics`： ```C++ defineDynamics( lattice, lattice.getBoundingBox(), new CylinderShapeDomain2D(cx, cy, radius), new BounceBack ); ``` 第二个参数是`Box2D`类型，它可用于提高此函数调用的效率：反弹节点归属于此框与`domainfunctional`指定的域之间的交集上的单元格。 如果您的计算域被指定为一组简单域的并集，您可以在每个更简单的域上迭代地使用`defineDynamics`。如果域是简单域的交集或任何其他几何运算，则需要在域泛函的定义中使用布尔运算符。 #### 布尔掩码的反弹域 在本节中，假设域的几何形状是从外部源规定的。 例如，当您模拟多孔介质并拥有有关多孔材料的扫描数据时，就会出现这种情况。这种最简单的方法包括将此数据存储为ASCII文件，该文件通过零和一区分流体节点和实体节点，并用空格、制表符或换行符分隔。该文件表示常规数组的内容，仅存储原始数据，没有关于矩阵维度的信息。数据布局必须与Palabos中的相同：z坐标的增量表示文件（在3D中）的连续进度，然后是y坐标，然后是x坐标。这些数据只是从文件中刷新到标量字段中。以下代码取自示例`examples/codesByTopic/io/loadGeometry.cpp`： ```C++ MultiScalarField2D boolMask(nx, ny); plb_ifstream ifile("geometry.dat"); ifile >> boolMask; ``` 请注意，目前没有为此I/O操作实施错误检查。您必须确保标量字段的维度与几何文件中的数据大小相对应。下一步，使用函数`defineDynamics`的其中一个版本实例化反弹节点： ```C++ defineDynamics(lattice, boolMask, new BounceBack, true); ``` 这个函数目前只为使用相同数据类型（在这种情况下为 T）的多块定义，这解释了为什么标量字段 `boolMask` 是T类型而不是人们所期望的bool类型。正如人们所期望的那样，标量场不是bool类型，而是真正的类型T。 ### 来自STL文件的反弹域 解析地定义计算域的形状并不总是很方便，甚至不可能通过任何现实的方式。在一般情况下，必须向Palabos提供描述计算域的表面网格。如果您决定通过BounceBack节点来描述域边界，如前几章所述，Palabos只会对域进行体素化，即将域的表面描述转换为体积描述，决定哪个流体节点（或体素）在域内。这种体素化是在Palabos内部完成的，并且是完全并行的：与整体模拟时间相比，这一步通常花费的时间可以忽略不计。 Palabos体素化器不仅决定哪些体素是流动的，哪些是反弹的，它还实施了一个内存节省方案，其中内存仅分配在流体覆盖的区域。 表面网格必须以二进制或ASCII STL格式提供给Palabos。许多表面网格格式可以使用开源程序`Meshlab`转换为STL。 可以在文件`examples/showCases/aneurysm/aneurysm_bounceback.cpp`中找到读取STL文件并创建由反弹节点包围的计算域的示例。 ### 来自STL文件的`off-lattice`边界条件 Palabos还提供了使用非晶格方案实现一般边界（入口、出口、壁面等）的可能性，这提供了比反弹节点的阶梯近似更好的精度。同样，域创建是并行且非常快的，并且体素器包括一种节省内存的域生成方法。 Palabos为BGK或MRT型流体以及热方程（对流-扩散）提供了一个全功能的非晶格框架。 目前还没有关于使用非晶格边界的完整文档。但是，可以在文件`examples/showCases/aneurysm/aneurysm.cpp`中找到一个示例。 ## 运行模拟 ### Palabos程序的时间周期 格子Boltzmann方法是一个显式求解器（或者说，Palabos中实现的最经典的格子Boltzmann方法）。一个迭代步后，流体的状态从$t$时刻过渡到$t+dt$时刻。接下来的讨论是在格子单位下进行的，并且$dt=1$。在Palabos中，格子Boltzmann方法的过程如下： 1. 开始时，所有的流体变量都是$t$时刻的，粒子分布函数处于“未碰撞”状态。 2. 在所有单元上进行碰撞步。粒子分布函数处于“碰撞后”状态。 3. 分布函数进行迁移。 4. 利用数据处理器处理非局部的或者不同单元之间耦合的数据。 5. 所有分布函数重新处于“未碰撞”状态，但是处于$t+1$时刻。 碰撞步（也即第2步），通过`lattice.collide()`方法实现，迁移步（也即第3步），通过`lattice.stream()`方法实现。这两种方法可以合并为`lattice.collideAndStream()`方法。在大多数平台上，`lattice.collideAndStream()`方法的计算效率更加高效，因为它只需要遍历计算域一次即可。 数据处理器可以通过调用`lattice.executeInternalProcessors()`方法手动调用，正如执行、集成和包装数据处理功能小节讲述的那样。然而，很少这样做，因为数据处理器会在函数`stream()`和函数`collideAndStream()`的后面自动执行。假如您需要Debug查找问题时，如果您不想调用完函数`stream()`和函数`collideAndStream()`后自动执行数据处理器，您可以使用这些函数的一些变体，如`lattice.collideAndStream(lattice.getBoundingBox())`。 在Palabos中，数据处理器总是在一个碰撞-迁移循环完毕后执行。因此，非本地操作总是在碰撞-迁移循环之后执行，此时分布函数处于“未碰撞”状态。但是，这不失一般性，因为在时间$t$的碰撞步之前执行操作相当于在时间$t-1$的迁移步之后执行它，对吗？唯一的问题出现在初始条件上，因为通常希望数据处理器在一开始就执行一次，就在设置初始条件之后。这是通过在开始第一个迭代步骤之前调用`lattice.initialize()`方法来实现的。该方法执行一次数据处理器，并在块内部执行一个内部通信步骤，以保证其内部状态一致。 ### 你在什么时候评估数据？ 为了监控程序的演变，不时评估一些流体力学变量是有用的，例如平均能量： ```C++ pcout << computeAverageEnergy(lattice) << endl; ``` 推荐您在碰撞步之前计算流体力学变量。虽然这种区别对于守恒变量如密度和速度并不重要（它们在碰撞前和碰撞后状态下相等），但对于非守恒变量（例如应力张量）很重要。非守恒速度矩只有在“未碰撞”状态下计算时才能与流体动力学变量相关。请注意，如果您使用`collideAndStream()`方法，则不会有做错事的风险，因为无论如何您都无法访问碰撞后变量。 通常，在调用方法`collideAndStream()`之后立即执行上例中的平均能量等流体力学变量的计算，因为此时所有流体力学变量都已明确定义，并且对应于同一时刻（碰撞和迁移过程中速度是明确定义的，但应变率不是）。格子的内部统计有一个例外。作为执行碰撞步骤的副作用，内部统计数据会自动计算而不会对性能产生任何影响（至少在串行程序中不会；对于并行情况，请参阅控制效率小节中的讨论）。然而，它们是在时间$t$时刻计算流体变量的，在碰撞-迁移循环期间，该循环将系统从时间$t$带到时间$t+1$。因此，通常在调用`collideAndStream()`方法后访问内部统计信息。在碰撞和迁移处理之前计算上例中的平均能量与从内部统计中访问平均能量产生相同的结果，如下面的示例所示，在碰撞和迁移处理之后： ```C++ pcout << getStoredAverageEnergy(lattice) << endl; ``` ### 其他需要做的重要事情 数字通常难以解释。因此，定期在程序中生成图像很有用，这样您就可以监控模拟状态，尽快发现问题，并在需要时重新运行程序。在二维和三维模拟中生成图像小节解释了如何执行此操作。 最后，不时保存模拟的状态总是好的，以免在计算机崩溃时丢失所有内容，并且如果您忘记生成重要的输出文件以供后期使用时能够恢复数据。这在检查点：保存和加载模拟状态小节中进行了解释。 ## 输入/输出 ### 输出流：向终端和文件中写入 在Palabos中，请永远不要使用`cout`,`cerr`,`clog`,`printf`等命令向终端输出。这些指令在Palabos中被其对应的并行版本`pcout`,`pcerr`,`pclog`代替（注意没有`pprintf`）。它们具有与对应的传统命令相同的语法和相同的行为，但此外，它们在并行程序中提供了预期的行为：如果您从在一百个内核上并行化的程序中打印消息，则该消息仅打印一次，而不是一百次。您可以使用它们将字符串、数字，甚至是格子、标量场或张量场的提取物打印到屏幕上： ```C++ pcout << "The average energy is " << computeAverageEnergy(lattice) << endl; pcout << "And here are some values of the energy, with 10-digit accuracy:" << endl; pcout << setprecision(10) << *computeEnergy(lattice, Box2D(0,10, 0,10)) << endl; ``` 同样的，数据可以流的方式传输到ASCII文件中（文件`examples/codesByTopic/io/useIOstream.cpp`中还提供了使用输出文件流的示例）。在这种情况下，标准`ofstream`被`plb_ofstream`替换，如下所示： ```C++ plb_ofstream ofile("energy.dat"); ofile << setprecision(10) << *computeEnergy(lattice) << endl; ``` 这对于随后将数据加载到交互式数据分析环境（如Octave或Matlab）中非常有用： ```C++ % Analysis of the data in Octave % Assign manually the dimensions: nx = 100; ny = 100; load energy.dat; energy = reshape(energy, nx, ny); imagesc(energy); ``` ### 输入流：从文件中读取大量数据 通过`plb_ifstream`类型的对象，输入文件的处理方式与上一节中的输出文件的处理方式相同。 例如，在前面的代码中写入的文件`energy.dat`可以通过以下命令流回矩阵： ```C++ plb_ifstream ifile("energy.dat"); if (ifile.is_open()) { // It is your responsibility to create the matrix with the right dimensions // nx-by-ny, compatible with the size of the data in "energy.dat". MultiScalarField2D energy(nx,ny); ifile >> energy; } ``` 要记住的重要一点是，必须首先使用正确的维度手动构建数据流入的标量字段，以避免内存错误。 然而，这种方法只有在数据被流式传输到Palabos对象时才有效。如果您希望将文件中的数据流式传输为纯`C++`数据类型，例如在访问用户定义的参数时，您可能不会使用如上所示的`plb_ifstream`数据类型，因为它会导致并行程序中出现意外行为。在这种情况下，存在两种可能的解决方法。第一种方法是推荐的方法，它要求用户以XML格式创建参数文件，然后使用Palabos的自动XML解析工具对其进行解析，如下面的从XML文件读取用户输入小节所述。 这样做有很多优点，因为它可以生成简短易读的代码，自动对用户输入进行类型转换，并在输入错误数据时自动进行错误处理。此外，XML格式的参数文件结构良好，并且在很大程度上是自动记录的。 如果由于某种原因您不愿意使用XML文件进行用户输入，您仍然可以使用普通的plb_ifstream文件，但您需要手动将数据传输到所有处理器以获得与数据并行兼容的工作并行程序Palabos的编程概念。为了说明这一点，让我们考虑一种情况，其中矩阵的维度nx和ny写在文件 energy.dat 的开头，并被读入程序以自动构造具有正确维度的标量场： ```C++ plb_ifstream ifile("energy.dat"); if (ifile.is_open()) { plint nx, ny; ifile >> nx >> ny; // Broadcast the input data to all processors; otherwise, the behavior // of the program is undefined. global::mpi().bCast(&nx, 1); global::mpi().bCast(&ny, 1); MultiScalarField2D energy(nx,ny); ifile >> energy; } ``` 代码`examples/codesByTopic/io/manualUserInput.cpp` 中说明了这种处理输入值的方式。请注意，与MPI相关的代码也在非并行程序中编译，在其中它根本不起作用。因此，最好在所有情况下都编写它们，以确保程序无需修改即可并行工作。 Palabos中没有用于交互式输入的对象pcin；必须通过输入文件或直接从命令行接收用户输入。 ### 获取命令行参数 一般在C++程序中，命令行参数可以通过主函数的`argc`和`argv`参数访问： ```C++ int main(int argc, char* argv[]) ``` 这在串行和并行程序中都可以正常工作。除此之外，Palabos还提供了两个名为`argc`和`argv`的全局对象，通过它们可以从程序中的任何位置访问命令行参数。一般建议使用这些全局对象而不是主参数，因为它们更智能，具有自动类型转换和错误处理。它们的用法在以下代码中说明： ```C++ int main(int argc, char* argv[]) { // Never forget to initialize Palabos, in every program. plbInit(&argc, &argv); pcout << "The number of input arguments is: " << global::argc() << endl; double double_argument; plint plint_argument; string string_argument; try { global::argv(1).read(double_argument); global::argv(2).read(plint_argument); global::argv(3).read(string_argument); } // React to errors, either if there are less than 3 arguments, // or if one of the arguments fails to convert to the expected // data type. catch(PlbIOException& exception) { // Print the corresponding error message. pcout << exception.what() << endl; // Terminate the program. exit(1); } // ... Execute the main program ... } ``` 在此示例中，C++异常机制用于对用户输入中的错误做出反应。如果您不想处理异常，可以使用`read`函数的`noThrow`版本： ```C++ bool ok = global::argv(1).readNoThrow(double_argument); if (!ok) { pcout << "Error while reading argument 1." << endl; } ``` 作为使用命令行参数的示例，您还可以查看程序`examples/showCases/rayleighTaylor2D.cpp`和`rayleighTaylor3D.cpp`。 ### 从XML文件读取用户输入 Palabos提供了一种从XML文件读取结构化输入数据的方法。此功能适用于串行和并行程序。然而，数据存储在普通的、非并行的数据结构中，这些数据结构在并行程序的所有线程上重复。因此，XML文件应仅用于小数据集，例如设置模拟所需的参数。大数据集，例如格子的内容，应该读入Palabos的数据结构，如输入流：从文件中读取大量数据小节所述，从而数据分布在并行机上。 XML数据在开源库`TinyXML`的帮助下进行解析。该库不需要手动安装，因为它与Palabos版本打包在一起。 Palabos可以读取任何格式良好的XML文档，但有以下三个限制： 1. 无法解析文档类型定义(DTD)或可扩展样式表语言(XSL)。 2. 无法识别属性（但即使存在属性，仍会正确解析标记）。例如，以下XML选项卡中的属性id值无法在Palabos中访问：``。 3. 给定的标签名称只能使用一次。如果多次使用，则在Palabos中只会访问最后一次出现的情况。 这些限制是为了简化Palabos中XML解析器的语法，并且因为它们不限制输入格式的通用性。以下是可与Palabos一起使用的典型输入文件： myInput.xml文件 ```xml /home/user/data/inputFile.dat 10 20 30 0.023 0.01 Dirichlet 0 0.1 0.2 0.3 0.4 0.4 0.3 0.2 0.1 0 ``` 以下是它们在Palabos程序中的解析方式： ```C++ try { // Open the XMl file. XMLreader xmlFile("myInput.xml"); // It's good policy to flush the content of the XML // file to the terminal, so that in future, when // you check the program's output, you remember // which parameters where used. pcout << "CONFIGURATION" << endl; pcout << "=============" << endl << endl; xmlFile.print(0); string inputFile; xmlFile["Geometry"]["inputFile"].read(inputFile); plint nx, ny, nz; xmlFile["Geometry"]["size"]["nx"].read(nx); xmlFile["Geometry"]["size"]["ny"].read(ny); xmlFile["Geometry"]["size"]["nz"].read(nz); T viscosity, uMax; xmlFile["Geometry"]["viscosity"].read(viscosity); xmlFile["Geometry"]["umax"].read(uMax); string inletKind; xmlFile["Inlet"]["kind"].read(inletKind); vector inletValues; xmlFile["Inlet"]["values"].read(inletValues); } catch (PlbIOException& exception) { pcout << exception.what() << endl; return -1; } ``` 特别要注意的是，可以从单个标签读取完整的数据数组，如上例中的最后一个标签所示，它被读入`inletValues`向量。 ### 二维和三维模拟中生成图像 可以从整个域、子域或标量场的切片生成图像。默认情况下，Palabos以PPM格式写入图像，这是一种ASCII格式，无需外部图形库即可轻松生成。如果安装了`ImageMagick`包（特别是这个包中包含的命令行工具`convert`），还可以以更标准的GIF格式写入图像。这在Linux和Windows下都有效。目录`examples/showCases`中的大多数示例说明了如何编写GIF图像。 作为第一步，您需要创建一个ImageWriter类型的对象，如下例所示： ```C++ ImageWriter imageWriter("leeloo"); ``` 参数字符串代表用于将标量变量映射到三分量RGB颜色方案的颜色图。可用的默认地图是地球、水、空气、火和leeloo。 下一步是编写PPM或GIF图像，要么使用调整到固定颜色值范围的颜色图，要么使用自动调整以适应数据的最小值和最大值的缩放范围。使用GIF格式，图像还可以重新缩放为不同的大小： ```C++ // Write a PPM image with a color map adapted to the range [minValue, maxValue]. imageWriter.writePpm(scalarField, "myImage", minValue, maxValue); // Write a PPM image with an automatically scaled color map. imageWriter.writeScaledPpm(scalarField, "myImage"); // Write a GIF image with a color map adapted to the range [minValue, maxValue]. imageWriter.writeGif(scalarField, "myImage", minValue, maxValue); // Write a GIF image with an automatically scaled color map. imageWriter.writeScaledGif(scalarField, "myImage"); // Write a GIF image with an automatically scaled color map, and rescale it in // order to fit in a siyeX-by-sizeY box (the aspect ratio is preserved, though). imageWriter.writeScaledGif(scalarField, "myImage", sizeX, sizeY); ``` 使用这种方法从三维模拟中生成快照也非常有用。在这种情况下，您需要从三维数据中提取一个切片，以具有一个退化维度（具有一个单元宽度的维度）的三维框的形式： ```C++ Box3D slice(0,nx-1, 0,ny-1, zValue, zValue); ImageWriter("earth").writeScaledGif(*computeVelocity(lattice, slice)); ``` ### 输出VTK文件进行后处理 来自块、标量场或张量场的所有数据都可以写入VTK数据格式的文件中。从那里，它可以通过科学数据可视化工具（如`Paraview`）进一步进行后处理。 Palabos中使用的VTK格式基于数据的无损二进制表示，通过`Base64`格式（它是一种基于ASCII的二进制格式，可能与您的电子邮件程序中的图像编码的格式相同）。使用保留模拟数据完整数值精度的格式通常是一种矫枉过正的做法，因为后处理操作通常需要的数值精度低于CFD模拟。为了节省一些内存，您至少可以将数据从双精度做单精度算术转换，如以下三维示例所示： ```C++ // Evaluate the discretization parameters, in order to write the data // in dimensionless units into the VTK file. T dx = parameters.getDeltaX(); T dt = parameters.getDeltaT(); // Create a VTK data object and indicate the cell width through the // parameter dx. This object is of the same type as the simulation, type T. VtkImageOutput3D vtkOut("simulationData.dat", dx); // Add a 3D tensor-field to the VTK file, representing the velocity, and rescale // with the units of a velocity, dx/dt. Explicitly convert the data to single- // precision floats in order to save storage space. vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt); // Add another 3D tensor-field for the vorticity, again as floats. vtkOut.writeData<3,float>(*computeVorticity(*computeVelocity(lattice)), "vorticity", 1./dt); // To end-with add a scalar-field for the density. vtkOut.writeData<1,float>(*computeDensity(lattice), "density", 1.); ``` 您可以向VTK文件添加任意数量的块，并且每个块都可以具有任意数量的组件（在上面的示例中，使用了3分量张量场和1分量标量场）。 类似地，可以使用`VtkImageOutput2D`类型的对象将二维数据写入VTK格式（示例在目录`examples/showCases/poiseuille/`中提供。然而，像Paraview这样的工具并不总是分析二维数据的最佳选择。通常建议您将数据刷新到文本文件中，如输出流：向终端和文件中写入小节中所述，然后使用Octave或Matlab（或 Python、R、SciLab或OO-Spreadsheet，或任何您喜欢的数据解释器将其可视化）。 ### 检查点：保存和加载模拟状态 虽然可以借助`VtkImageOutputXD`对象将模拟结果写入VTK文件，但没有实现相反的过程：无法将VTK数据读入Palabos的数据结构。如果您希望保存系统状态以便稍后恢复模拟，请使用函数`saveBinaryBlock`，例如文件`examples/codesByTopic/io/checkPointing.cpp` 中的说明： ```C++ saveBinaryBlock(lattice, "checkpoint.dat"); ``` 并稍后使用函数`loadBinaryBlock`加载它： ```C++ loadBinaryBlock(lattice, "checkpoint.dat"); ``` 使用此过程，您一次只能写入一个数据块。如果模拟的状态由各种块表示（例如在使用几个块格的多组分流体中），则每个块都需要保存在单独的文件中。目前无法保存结构数据（动态对象和数据处理器）。因此，在加载模拟状态时，必须在加载数据之前重新创建上下文，例如边界条件。 ## 数据评估 ### 概述 在格子Boltzmann方法中，我们模拟的变量主要是分布函数。这与传统的CFD工程师所关注的宏观变量，如密度、压力、速度等形成了鲜明的对比。显然，在将数据提供给标准的后处理工具之前，需要以某种方式转换数据。Palabos中提供了许多用于转换、评估或转换数据的函数，如附录中列出的用于数据分析和其他目的的非可变操作。在本节中，以函数`computeVelocity`为例进行讨论，因为其他函数的工作方式非常相似。这些函数是为二维和三维情况定义的，它们适用于单子块和多块。为了便于说明，以下代码参考了三维多块的情况。 函数`computeVelocity`有三个版本： ```C++ // Version 1 void computeVelocity(MultiBlockLattice3D& lattice, MultiTensorField3D::d>& velocity, Box3D domain); // Version 2 std::unique_ptr > computeVelocity(MultiBlockLattice3D& lattice, Box3D domain); // Version 3 std::unique_ptr > computeVelocity(MultiBlockLattice3D& lattice); ``` 在第一个版本中，速度是在块格的子域上计算的，并将结果写入三分量张量场的相应子域。块格和张量场不需要具有相同的大小，也不需要具有相同的内部块排列。 如果域超过块格或张量场的维度，则相应地修剪域。 在实践中更常用的第二个版本中，自动创建了一个域大小的张量场，并从函数中返回。函数返回值中使用的`auto-pointer`关键字是指`C++`标准库提供的一类智能指针。从用户的角度来看，它实际上与指针相同。 这意味着您对待`std::unique_ptr >`类型的对象的方式与对待`MultiTensorField3D*`类型的对象相同。两者最大的区别在于自动指针有一个自动内存管理机制，你永远不需要对这种类型的指针调用操作符delete。如果函数的返回值是一个原始指针，您将无法对不同的运算符进行流水线操作，如下一节所示，因为您永远无法获得指向它们的显式指针，因此无法删除它们。 以下代码显示了函数`omputeVelocity`的典型用例： ```C++ pcout << *computeVelocity(lattice, domain) << endl; ``` 在这里，计算出的速度值会立即重定向到终端。函数调用前面的星号用于取消对速度场的指针的引用。在此程序行的末尾，速度场的内存会自动释放，不需要显式释放。 第三个版本是一个纯便利函数，其中参数域被格子的完整域`lattice.getBoundingBox()`替换。 ### 流水线数据评估运算符 像`computeVelocity`这样的函数的返回值可以直接用作其他数据评估运算符的参数。这样，就可以构造复杂的表达式。例如，上一节中的函数`computeAverage`可以替换为计算速度场，然后计算每个元素的范数平方，除以二，最后计算平均值： ```C++ pcout << "The value " << *computeAverageEnergy(lattice) << " is the same as " << computeAverage ( *multiply ( 0.5, *computeNormSqr (*computeVelocity(lattice) ) ) ) << endl; ``` 本示例中使用的所有函数都列在附录附录B 部分功能/类的参考中。目录`examples/codesByTopic/scalarField`中提供了更多数据评估、标量字段构造和数据评估运算符组合的示例。 ## 粒子 ### 概览 Palabos为模拟粒子和实现粒子-流体相互作用提供了一个完全并行的框架。我们目前正在处理此功能的文档。 同时，您可以通过查看示例代码`codesByTopic/particlesInTube`获得第一印象。 ## 网格加密 ### 概览 从0.7版本之后，Palabos支持网格加密功能。目前，只实现了二维问题的网格加密功能。至于Palabos的其他方面，网格加密与其他实现的概念是正交的。这意味着，例如，以前实现的动态对象和数据处理器也可以与网格加密一起使用。此外，加密网格可以以一般方式并行化（可以不受限制地选择不同进程上的域之间的边界）。 ### 多层网格加密 Palabos中的网格加密是多块概念的扩展，被称为“多网格”。每一级网格加密都由一定数量的块覆盖，并由多块表示。完整的网格——“多网格”——由所有这些层次的叠加构成。 以二维顶盖驱动腔为例，由三层细网格表示，以提高顶盖附近的精度，特别是左上角和右上角附近的精度： ![alt](./figs/multi-block.jpg) 诀窍是：Palabos只需为您生成三个单独的多块，每个块都将数据保存在一定程度的网格加密。多块具有稀疏的内存表示，这意味着整个操作在内存和效率方面不会花费任何成本。在我们的例子中，有三个多块，如下图所示。最左边的多块是最小的一次，因为它保存了最粗级别的数据，最右边的最大的代表最精细的级别。 ![alt](./figs/multi-block2.jpg) 然后，Palabos将三个多块聚集到一个数据结构中（但为方便起见，您仍然可以单独访问它们），称为“多网格”结构： ![alt](./figs/multi-block3.jpg) #### 交互式创建多网格 示例程序`showCases/gridRefinement2d`说明了多重网格数据结构的用法。 ## 并行编程 ### 并行编程方法 使用Palabos编写的程序可以自动并行化，至少如果它们遵循本用户指南中提供的建议的话。并行化是使用MPI库的消息传递范式执行的。这种方法适用于分布式内存平台（例如集群）以及共享内存平台（SMP或多核架构）。事实上，由于双核或多核处理器的高可用性，建议几乎总是并行编译Palabos程序，即使在开发阶段也是如此。通过利用机器的所有内核获得的显着加速可以提高代码开发和程序测试周期的效率。 为了实现具有串行应用程序外观和感觉的程序的并行性，Palabos区分了两类数据。空间分布的数据，例如块格、标量或张量场，通过数据并行范式处理。数据空间被划分为较小的区域，这些区域分布在并行机器的节点或核心上。其他需要少量存储空间的数据类型在每个节点上都有复制。借助 MPI 的单程序多数据模型，所有本机C++数据类型以及C++ STL的数据类型和Palabos的小数据类型（如`Array`）都会自动复制。这些类型应该用于处理标量值，例如模拟的参数或解决方案的积分值（例如平均能量）。 最后，再次指出，只有多块结构（多块点阵、多标量场和多张量场）是数据并行的，而它们的原本的对应物如果在终端中创建，则会被复制到所有核心中。出于这个原因，为了获得并行性能改进，始终使用多块结构非常重要。 ### 控制效率 #### 高效并行代码遵守的规则 用户指南中多次强调并重复在Palabos中高效程序实施的首要原则，因为它非常重要。 --- **高效编程的首要原则**： 最终用户应用程序永远不应包含在多块索引上运行的循环。偶尔可以通过类似数组的接口访问多块的单元格，以评估或修改数据。但是在空间循环中重复访问它们是不行的，因为这在并行程序中会非常缓慢。 空间循环必须出现在数据处理器内部而不是其他任何地方。 --- 一旦遵守了这条规则，导致高效并行程序的大多数其他要素都是简单的常识。例如，输入/输出操作在并行程序中相对较慢，因为它们当前的实现方式（参见瓶颈部分）。因此，您应该确保程序不会过于频繁地输出数据。例如，目录`examples`中提供的许多程序在并行时非常慢，因为它们经常生成图像来说明一个点。如果您尝试对它们的并行性能进行基准测试，那么您应该注释掉输出操作。 在大多数情况下，还可以通过避免集体MPI操作（例如归约）来显着提高程序的效率。格子Boltzmann方法实际上并不需要归约操作来执行它们的时间迭代。出于方便，在每个时间步都执行缩减，以评估内部统计数据（平均能量和平均密度）。这些内部统计数据在串行程序中非常有用，因为它们实际上是免费计算的，但它们实际上会破坏并行程序的性能。这是因为减少就像一个同步屏障，这一事实被观察到对性能有负面影响，尤其是在类似集群的并行机器上。因此，可以关闭内部 统计信息（以及归约操作），通过函数调用 ```C++ lattice.toggleInternalStatistics(false); ``` 如果您只是想不时了解内部统计信息，例如，二十次迭代一次，可以将它们打开一小段时间， ```C++ lattice.toggleInternalStatistics(true); ``` 然后，执行时间迭代： ```C++ lattice.collideAndStream(); ``` 最后，访问统计信息并再次将其关闭 ```C++ T rho = getStoredAverageDensity(lattice); lattice.toggleInternalStatistics(true); ``` #### 使用多少个处理器？ 当一个程序在少数核上并行化时，并行效率通常接近最优。内核数量每增加一倍，执行速度也几乎翻倍。随着内核数量的进一步增加，加速曲线开始变平，并且在某些时候不值得添加更多内核。通常将程序的效率称为一个内核上的执行时间，除以给定数量n个内核上的执行时间，再除以n。如果程序以75%的效率运行，你可以说它运行得相当好，因为并行化从来都不是最优的。另一方面，如果效率下降到50%以下，则要求您减少内核数量是合理的，因为您只是在浪费研究所的资源。 在不深入讨论并行理论的细节的情况下，让我们只提一下用于解决问题的最佳内核数量取决于许多参数。对于大型域的模拟，并行效率往往更好，并且受到良好互连网络的有利影响（实际上是通信速度与核心计算速度的高比率）。在实践中，必须根据经验确定最佳核心数 对于给定的问题和给定的并行机器。一个容易记住的值是当程序在最佳状态下运行时与每个内核相关联的子域的大小。例如，在给定的本土集群上，当每个内核域大小（总域大小除以内核数）约为 20×20×20 时，观察到示例程序`cavity3d`以75%的效率运行。从那时起，这台机器与“20立方规则”相关联，通过它可以轻松快速地评估用于给定问题规模的最佳核心数量。 #### 当前实现中的瓶颈 尽管投入大量精力以在Palabos中获得出色的并行效率，但仍有改进的余地。目前，最严重的瓶颈之一似乎是输入/输出操作。对于输入和输出，数据始终通过其中一个核心传输。因此，这些操作是非并行的：它们并不快，或者在多核上比在单核上快一点。虽然这是可以的 在具有数百个内核的较小并行机器上运行，它往往会在千核硬件上出现问题，因为输入/输出操作可以支配程序的整体执行时间。因此，并行输入/输出是需要在不久的将来实现的最重要的功能之一。 ## 用于非本地操作和分块耦合的数据处理器 动力学类和数据处理器是Palabos程序的基本要素。从概念上讲，您可以将动力学类视为原始“格子Boltzmann范式”的实现，而数据处理器则是更通用的数据并行多块范式的实现。动力学类（相对）容易，因为格子Boltzmann方法很容易，只要没有边界、细化网格、并行程序或任何其他高级结构成分。数据处理器可能有点难以理解，因为它们有一项艰巨的任务是解决动力学类不负责的“其他所有问题”。它们实现了模型的所有非局部成分，它们在标量场和张量场上执行操作，并在所有类型的块之间创建耦合。而且，就像动态对象一样，它们处理归约操作，它们必须高效且本质上可并行化。 当您使用标量场和张量场时，数据处理器的重要性尤其明显。与块格不同，这些数据字段没有智能单元，其具有对动力学对象的本地引用；因此，数据处理器是对其执行集体操作的唯一有效方式。举个例子，考虑一个表示速度场的`TensorField3D` 类型的场（每个元素都是`Array`类型的速度向量）。例如，通过调用函数`computeVelocity` 从LB模拟中获得这样的字段。通过有限差分计算速度的空间导数可能很有趣。这样做的唯一道德上正确的方法是通过数据处理器，因为它是唯一可并行化、可扩展、高效且与未来版本的`Palabos`前向兼容的方法。这正是Palabos在调用`computeVorticity()` 或`computeStrainRate()`等函数之一时所做的事情（参见附录张量场->张量场）。请再次注意，您还可以通过在张量场的空间索引上编写一个简单的循环来评估有限差分方案。这将在串行和并行中产生正确的结果，并且在串行中甚至会非常有效，但它不是可并行的（并行时，效率是损失而不是获得）。 总而言之，应该清楚的是，虽然数据处理器是强大的对象，但它们也需要解决一定程度的复杂性，因此在概念上不如Palabos的其他组件简单。因此，数据处理器肯定是Palabos用户界面中最难理解的部分，而这一部分充满了临时规则，您只需要在某个时候学习。为了稍微缓解这个问题，本节从两个相对简单的主题开始，它们教你如何通过可用的辅助函数在Palabos中解决许多任务，并避免在许多情况下显式编写数据处理器。本节的其余部分包含许多指向Palabos中现有实现的链接，强烈建议您实际查看这些示例以吸收理论概念。 ### 使用辅助函数来避免显式编写数据处理器 根据使用它们的用途，数据处理器可以（任意）分为三类。第一类是关于设置模拟、将动力学对象分配给格子的子域、初始化分布函数等。这些方法在初始化速度和密度和设定边界条件部分中进行了解释，并且这些函数在附录用于模拟设置和其他目的的可变（就地）操作小节有介绍。第二类包括添加到格子中以实现物理模型的数据处理器。许多模型是预定义的，如实现的流体模型小节所述。最后一类是用于评估数据，并用于执行简短的后处理任务，如上面计算速度梯度的示例。示例在数据评估小节中找到，可用函数列在附录用于数据分析和其他目的的非可变操作中。 ### 一些对于本地操作很便利的包装器 想象一下，您必须在一个块格子上执行本地初始化任务，即您不需要访问周围邻居的值的任务，并且可用的Palabos函数不足。作为编写经典数据处理函数的轻量级替代方案，您可以实现一个继承自`OneCellFunctionalXD`并实现虚拟方法的类（此处为二维问题） ```C++ virtual void execute(Cell& cell) const; ``` 在此方法的主体中，只需对参数单元格执行所需的操作。如果动作取决于空间位置，则可以改为从`OneCellIndexedFunctionalXD`继承并实现该方法（再次针对二维情况进行说明） ```C++ virtual void execute(plint iX, plint iY, Cell& cell) const; ``` 然后通过函数调用将这个函数的一个实例应用于格子，例如 ```C++ // Case of a plain one-cell functional. applyIndexed(lattice, domain, new MyOneCellFunctional); // Case of an indexed one-cell functional. applyIndexed(lattice, domain, new MyOneCellIndexedFunctional); ``` 此方法用于位于`examples/showCases/multiComponent2d`的示例程序中。在这里，需要一个定制的初始化过程来访问格子的外部标量以进行热模拟。 这些仅用于一个格子的函数不如通常的数据处理泛函通用，因为它们不能用于非本地操作。此外，它们在数值上往往效率较低，因为Palabos需要对域的每个单元上的单元功能的方法执行虚函数调用。然而，这种效率损失通常在初始化阶段完全可以忽略不计，在这个阶段，拥有可在并行机器上扩展的代码很重要，而不是优化代码以获得微秒的增益。在这种情况下，您应该更喜欢仅用于一个格子的函数所建立的更短且更易于阅读的代码。 ### 编写数据处理器（或者说，是数据处理函数） 对矩阵执行操作的一种常见方法是在矩阵的所有元素上或至少在给定的子域上编写某种循环，并对每个单元格执行操作。如果矩阵的内存被细分为更小的组件，就像Palabos的多块的情况一样，并且这些组件分布在并行机器的节点上，那么您的循环也需要细分为相应的更小的循环。Palabos中的数据处理器的目的是为您自动执行此细分。然后它提供细分域的坐标，并要求您在这些域上执行操作，而不是原始的完整域。 作为数据处理器的开发人员，您几乎总是与所谓的数据处理功能保持联系，这些功能通过在幕后执行部分重复性任务来提供简化的界面。存在许多不同类型的数据处理函数，如下一节所述。为了说明起见，我们现在考虑一个数据处理器的情况，它作用于单个二维单块点阵或多块点阵，并且不执行任何数据缩减（它不返回任何值）。假设该操作的目的是在格子单元上交换给定数量的`f[1]`和`f[5]`的值。这可以（并且应该，为了代码清晰起见）通过一个简单的仅用于一个格子的函数来完成，该函数在一些对于本地操作很便利的包装器小节中介绍，但我们现在想做真正的事情，并进入数据功能的核心。 该作业的数据处理函数必须继承自`BoxProcessingFunctional2D_L`（L 表示数据处理器作用于单个格子），并实现以下描述的其他方法，虚拟方法过程： ```C++ template class Descriptor> class Invert_1_5_Functional2D : public BoxProcessingFunctional2D_L { public: // ... implement other important methods here. void process(Box2D domain, BlockLattice2D& lattice) { for (plint iX=domain.x0; iX<=domain.x1; ++iX) { for (plint iY=domain.y0; iY<=domain.y1; ++iY) { Cell& cell = lattice.get(iX,iY); // Use the function swap from the C++ STL to swap the two values. std::swap(cell[1], cell[5]); } } } }; ``` 函数过程的第一个参数对应于由Palabos在细分原始域以适应原始块的小组件后计算的域。第二个参数对应于要在其上执行操作的格子。这个参数总是一个原子块（即它总是一个块格，而不是多块格），因为在函数调用过程的时候，Palabos已经细分了原始块并正在访问它的内部。如果您将此与编写一个单元功能的过程进行比较，如一些对于本地操作很便利的包装器小节所示，您会发现在当前情况下您需要做的额外工作是自己编写空间索引上的循环的域。一直手动写出这些循环很累，尤其是当您编写许多数据处理器时，而且容易出错。但这是为获得最佳效率而付出的代价，在计算物理领域，效率比软件工程的其他领域更重要，不幸的是，有时必须重视效率，而不是优雅。 与仅用于一个格子的功能相比的另一个优势是可以实现非本地操作。在以下示例中，`f[1]`的值与右侧相邻单元格上的`f[5]`交换： ```C++ void process(Box2D domain, BlockLattice2D& lattice) { for (plint iX=domain.x0; iX<=domain.x1; ++iX) { for (plint iY=domain.y0; iY<=domain.y1; ++iY) { Cell& cell = lattice.get(iX,iY); Cell& partner = lattice.get(iX+1,iY); std::swap(cell[1], partner[5]); } } } ``` 如果您遵守两条规则，您可以这样做，而不会冒险访问格子范围之外的单元格： 1. 在最近邻格子（D2Q9、D3Q19 等）上，您可以是一个单元格的非本地格子，但不能更多（您可以写`lattice.get(iX+1,iY)` 但不能写`lattice.get(iX+2, iY)`。在具有扩展邻域的格子上，您还可以扩展访问数据处理器中相邻单元的距离。允许的非局部性数量由常量`Descriptor::vicinity`确定。 2. 作用于通信包络的数据处理器中不允许进行非本地操作（您需要覆盖的方法小节解释了这意味着什么）。 为了结束本小节，让我们总结一下您可以在数据处理器中做的事情，以及不允许做的事情。您可以访问提供的范围内的单元格（加上最近的邻居，或根据格子拓扑结构的更多邻居），读取它们并修改它们。您执行的操作可以依赖于空间，但这种空间依赖必须是通用的，不能依赖于提供给函数进程的参数域的特定坐标。这是非常重要的一点，我们将其列为数据处理器的首要规则： --- **数据处理器的首要规则**： 数据处理器必须始终以这样一种方式编写，即在给定域上执行数据处理器与将域分成两个子域，然后在每个子域上连续执行数据处理器具有相同的效果。 --- 实际上，这意味着您不能根据参数域的参数`x0`、`x1`、`y0`或`y1`或直接根据索引`iX`或`iY`做出任何逻辑决策。相反，这些局部索引必须首先转换为全局索引，独立于数据处理器的细分，如绝对位置和相对位置小节所述。 #### 数据处理功能的目录 根据应用数据处理器的块的类型，存在不同类型的处理功能，如下所示： **Class BoxProcessingFunctionalXD** void processGenericBlocks(BoxXD domain, std::vector*> atomicBlocks); 这个类实际上从未使用过。当其他一切都失败时，它是备用选项。它处理任意数量的任意类型的块，这些块被强制转换为通用类型`AtomicBlockXD`。在使用之前，您需要将它们转换回它们的真实类型。 **Class LatticeBoxProcessingFunctionalXD** void process(BoxXD domain, std::vector*> lattices); 使用此类来处理任意数量的块格，并可能在它们之间创建耦合。例如，该数据处理函数用于为在文件`src/multiPhysics/shanChenProcessorsXD.h`和`.hh`中定义的Shan/Chen多分量模型定义任意数量的格子之间的耦合。这种类型的数据处理函数也不是很常用，因为下面列出的两个版本在大多数情况下更合适。 **Class ScalarFieldBoxProcessingFunctionalXD** void process(BoxXD domain, std::vector*> fields); 和上面一样，用于标量场。 **Class TensorFieldBoxProcessingFunctionalXD** void process(BoxXD domain, std::vector*> fields); 和上面一样，用于张量场。 **Class BoxProcessingFunctionalXD_L** void process(BoxXD domain, BlockLatticeXD& lattice); 作用于单个格子的数据处理器。 **Class BoxProcessingFunctionalXD_S** void process(BoxXD domain, ScalarFieldXD& field); 作用于单个标量的数据处理器。 **Class BoxProcessingFunctionalXD_T** void process(BoxXD domain, TensorFieldXD& field); 作用于单个张量的数据处理器。 **Class BoxProcessingFunctionalXD_LL** void process(BoxXD domain, BlockLatticeXD& lattice1, BlockLatticeXD& lattice2); 用于处理耦合具有潜在不同描述符的两个格子的数据处理器。类似地，对于两个标量场有一个`SS`版本，对于两个具有潜在不同维度nDim的张量场有一个`TT`版本。 **Class BoxProcessingFunctionalXD_LS** void process(BoxXD domain, BlockLatticeXD& lattice, ScalarFieldXD& field); 用于处理耦合格子和标量场的数据处理器。类似地，对于格子张量和标量张量情况，有一个`LT`和一个`ST`版本。 对于这些处理函数中的每一个，对于数据处理器执行归约操作并返回值的情况，存在一个“归约”版本（例如`ReductBoxProcessingFuncionalXD_L`）。 #### 您需要覆盖的方法 除了方法`process`之外，数据处理功能必须重写三个方法。现在以本节开头介绍的类`Invert_1_5_Functional2D`为例说明这三种方法的使用： ```C++ BlockDomain::DomainT appliesTo() const { return BlockDomain::bulk; } void getModificationPattern(std::vector& isWritten) const { isWritten[0] = true; } Invert_1_5_Functional2D* clone() const { return new Invert_1_5_Functional2D(*this); } ``` 首先，您需要提供方法`clone()`，它在Palabos中是典型的（参见使用Palabos编写代码一节）。接下来，您需要告诉Palabos数据处理器处理的哪些块正在被修改。在本例中，只有一个块。在一般情况下，向量`isWritten`的大小等于所涉及块的数量，并且您必须为每个块分配一个标志`true/false`。其中，Palabos利用此信息来决定在数据处理器执行后是否需要块的进程间通信。 第三种方法`applysTo`方法用于决定数据处理器是仅作用于模拟的实际域（BlockDomain::bulk）还是还包括通信包络（BlockDomain::bulkAndEnvelope）。让我们记住，作为多块组件的原子块由单个单元层（或扩展格的多个单元层）扩展，以合并块之间的通信。这个包络与另一个原子块的大部分重叠，并且信息在相应的块和包络单元之间复制。正是这个信封使得实现非本地数据处理器成为可能，而不会招致访问超出范围数据的危险。通常，执行大量原子块上的数据处理器就足够了，最好这样做，以避免在使用信封时出现下面列出的限制。这已经足够了，因为信封和大量相邻块之间会自动启动通信，以便在需要时更新信封中的值。仅当 (1) 将新的动力学对象分配给部分或所有单元格时才需要包含包络（就像调用函数`defineDynamics`时所做的那样），或 (2) 如果您修改动力学对象的内部状态（就像调用函数`defineVelocity`以在`Dirichlet`边界节点上分配新的速度值时所做的那样）。在这些情况下，包括包络是必要的，因为动力学对象的性质和内容在原子块之间的通信步骤中不会被传输。传输的唯一信息是单元数据（分布函数的数量和外部标量）。 如果您决定将信封包含在数据处理器的应用程序区域中，则必须遵守以下两个规则。否则，将出现未定义的行为。 1. 数据处理器必须完全是本地的，因为没有额外的信封可用于处理非本地数据访问。 2. 数据处理器最多可以对所涉及的块之一具有写访问权（从方法`getModificationPattern()`返回的向量`isWritten`最多可以在一个位置具有值`true`）。 #### 绝对位置和相对位置 在数据处理器的空间循环中使用的坐标`iX`和`iY`除了执行循环外，对其他任何事情都毫无用处，因为它们代表原子块的局部变量，原子块本身位于整个多块内部的随机位置。因此，要根据空间位置做出决定，必须首先将局部坐标转换为全局坐标： ```C++ // Access the position of the atomic-block inside the multi-block. Dot2D relativePosition = lattice.getLocation(); // Convert local coordinates to global ones. plint globalX = iX + relativePosition.x; plint globaly = iY + relativePosition.x; ``` 目录`examples/showCases/boussinesqThermal2d/`中提供了一个示例。这种转换有点尴尬，这也是使用仅作用于一个格子的函数的一个很好的理由，该函数在一些对于本地操作很便利的包装器小节中介绍，它会自动为您完成这项工作。 类似地，如果您在多个块上执行数据处理器，则所有相关块中的相对坐标不一定相同。如果您在全局坐标中测量事物，那么方法`process`的参数域总是与所有涉及的重叠块。这是由在Palabos中实现的算法所保证的。然而，应用数据处理器的所有多块不一定使用相同的内部数据分布，并且可能对局部坐标有不同的解释。方法`process`的参数域总是作为第一个原子块的局部坐标提供。要获得其他块的坐标，必须应用相应的转换： ```C++ Dot2D offset_0_1 = computeRelativeDisplacement(lattice0, lattice1); Dot2D offset_0_2 = computeRelativeDisplacement(lattice0, lattice2); plint iX1 = iX + offset_0_1.x; plint iY1 = iY + offset_0_1.y; plint iX2 = iX + offset_0_2.x; plint iY2 = iY + offset_0_2.y; ``` 同样，此过程在`examples/showCases/boussinesqThermal2d/`中的示例中进行了说明。如果验证了以下任一条件，则需要计算此位移（如果不确定，最好默认计算位移）： 1. 应用数据处理器的多块不具有相同的数据分布，因为它们的构造不同。 2. 应用数据处理器的多块没有相同的数据分布，因为它们的大小不同。这适用于所有函数，例如`computeVelocity`，它计算格子域上的速度。它使用一个数据处理器，作用于原始格子（很大）和速度场（可以更小，因为它具有子域的大小）。 3. 数据处理器包括信封。在这种情况下，相对位移源于块节点与来自不同原子块的包络节点耦合的事实。这也是为什么通常最好不要将信封包含在数据处理器的应用程序域中的另一个原因。 #### 执行、集成和包装数据处理功能 基本上有两种使用数据处理器的方法。在第一种情况下，处理器通过调用函数`executeDataProcessor`仅在一个或多个块上执行一次。在第二种情况下，通过调用函数`addInternalProcessor`将处理器添加到块中，然后承担内部数据处理器的角色。内部数据处理器是块的一部分，可以通过调用该块的`executeInternalProcessors`方法多次执行。当数据处理步骤是流体求解器算法的一部分时，通常会选择这种方法。例如，考虑边界条件的非局部部分、多组分流体中的组分之间的耦合，或使用`Boussinesq`近似的热代码中流体与温度场之间的耦合。在块格中，内部处理器具有特殊的作用，因为方法`executeInternalProcessors`在方法`collideAndStream()`和方法`stream()`的末尾自动调用。这种行为是基于这样的假设：`collideAndStream()`代表一个完整的格子Boltzmann迭代循环，并且如果使用`stream()`，它位于该循环的末尾。因此，内部处理器被认为是格子Boltzmann迭代的一部分，并在碰撞和迁移步骤之后的最后执行。 为方便起见，在数据处理功能的目录小节介绍的每种数据处理函数重新定义了对`executeDataProcessor`和`addInternalProcessor`的函数调用，新函数分别称为`applyProcessingFunctional`和`integrateProcessingFunctional`。例如，要在给定格子和标量字段的整个域上执行`BoxProcessingFunctional2D_LS`类型的数据处理函数（它们可以是多块或原子块类型），要使用的函数调用具有以下形式： ```C++ applyProcessingFunctional ( new MyFunctional, lattice.getBoundingBox(), lattice, scalarField ); ``` Palabos中所有预定义的数据处理函数都额外包装在一个便利函数中，以简化语法。例如，用于二维字段的函数`computeVelocityNorm`的三个版本之一在文件`src/multiBlock/multiDataAnalysis2D.hh`中定义如下： ```C++ template class Descriptor> void computeVelocity( MultiBlockLattice2D& lattice, MultiTensorField2D::d>& velocity, Box2D domain ) { applyProcessingFunctional ( new BoxVelocityFunctional2D, domain, lattice, velocity ); } ``` #### 内部数据处理器的执行顺序 在函数调用`executeInternalProcessors()`中有不同的方法来控制内部数据处理器的执行顺序。首先，每个数据处理器都归属于一个处理器级别，这些处理器级别以递增顺序遍历，从级别0开始。默认情况下，所有内部处理器都归属于级别0，但您可以将它们放入任何其他级别，指定为函数`addInternalProcessor`或`integrationProcessingFunctional`的最后一个可选参数。在处理器级别内，数据处理器按照它们添加到块中的顺序执行。除了强制执行顺序之外，数据处理器对给定级别的归属对多块内的通信模式有影响。事实上，通信不是在执行具有写访问权限的数据处理器之后立即执行的，而是仅在从一个级别切换到下一个级别时执行。这样，一个级别内的数据处理器所需的所有MPI通信都被捆绑并更有效地执行。为了澄清这种情况，让我们写下块格的一个迭代周期的细节，该块格的数据处理器位于级别0和级别 1，并在函数调用`collideAndStream`结束时自动执行它们： 1. 执行本地碰撞，然后执行迁移步骤。 2. 在级别0执行数据处理器。目前尚未进行任何通信。因此，该级别的数据处理器执行非本地操作的能力有限，因为通信包络中的单元数据是错误的。 3. 在块格的原子块之间执行通信以更新包络。如果任何其他外部块（格子、标量场或张量场）被级别0的任何数据处理器修改，则也要更新这些块中的包络。 4. 在级别1执行数据处理器。 5. 如果块格或任何其他外部块被级别1的任何数据处理器修改，则相应地更新包络。 虽然这种行为可能看起来有点复杂，但它会导致程序的直观行为，并提供一种控制数据处理器执行的通用方法。需要特别强调的是，如果一个数据处理器B依赖于另一个数据处理器A先前产生的数据，则必须确保A和B之间存在适当的因果关系。在所有情况下，B必须在A之后执行。此外，如果B是非本地的（因此访问信封上的数据）并且A是仅批量数据处理器，则需要在A和B的执行。因此，必须在不同的处理器级别上定义A和B。 如果您手动执行数据处理器，您可以选择仅执行给定级别的处理器，方法是将级别指示为方法`executeInternalProcessors(plint level)` 的可选参数。还应该提到，处理器级别可以具有负值。负处理器级别的优点是它不会通过默认函数调用`executeInternalProcessors()`自动执行。只能通过调用`executeInternalProcessors(plint level)`手动执行。对于经常执行但不是在每个迭代步骤中执行的数据处理器，利用这种行为是有意义的。每次调用`applyProcessingFunctional`的效率会稍低一些，因为数据处理器在内部原子块上的分解会产生开销。 ## 实用程序 ### 计时器 Timer对象可用于进行时间测量和性能评估。它们根据挂钟时间测量时间间隔，在串行程序中借助`C`函数`clock()`，在并行程序中借助`MPI`函数`MPI_Wtime()`。通过类似的命令访问计时器对象 ```C++ global::timer("nameOfTimer") ``` 其中`nameOfTimer`是您选择的任意字符串。您可以根据需要使用尽可能多的不同计时器对象，并且在您第一次引用它们时会自动创建计时器。Timer对象是全局的，可以累积测量时间间隔。因此，您可以开始测量程序的一个文件中的时间间隔，然后通过引用具有相同字符串参数的计时器，从另一个文件中的某个位置累加更多时间。 计时器的行为类似于秒表。当您调用方法`global::timer("nameOfTimer").start()`时，它开始测量时间。调用`stop()`方法后，时钟停止移动，但会从下一次调用`start()`时的位置开始。要将时钟重置为零，请调用方法`reset()`或`restart()`。并且在目录`examples/codesByTopic/smagorinskyModel`中提供了使用计时器对象的示例。 ### 数据追踪 数据追踪器最常用于确定模拟何时达到稳定状态。他们跟踪标量值（例如平均能量）的时间演变，并确定当在固定时间窗口上测量的该量的标准偏差低于给定阈值时达到稳定状态。目录`examples/showCases/boussinesqThermal2D`和对应的3D目录中提供了一个示例。 ## 附录 ### 附录A 示例程序 目录`examples/tutorials`中的内容将在Palabos Tutorial中单独介绍。 请注意：这些示例代码的唯一目的是教您使用Palabos，在许多情况下，我们没有从科学的角度验证结果。如果您在研究中使用代码，您有责任理解并可能调整代码。 #### 目录`examples/showCases` 该目录包含用于模拟众所周知的基准问题的全功能（但基本）程序。 ##### aneurysm 这是功能最全面的示例，因此不是开始使用Palabos的最佳选择。它说明了高级Palabos模拟的关键概念： 1. 读取STL格式的网格，对域进行体素化，并生成off-lattice边界条件。 2. 将模拟结果从粗网格复制到细网格。在此示例中，该技术用于加速收敛到稳态（首先在粗网格上收敛，并将其用作细网格上的初始条件）。 3. 从XML文件中读取用户输入参数。 4. 各种类型的输出，包括表面数据（壁面剪应力）。 ##### boussinesqThermal2d 限制在热底壁（速度无滑移）和冷顶壁（速度无滑移）之间的流体。侧壁是周期性的。在重力的影响下，形成热对流。热效应通过Boussinesq近似建模：流体是不可压缩的，温度的影响只能通过体力项可见，代表浮力效应。温度场服从对流-扩散方程。 首先以完全对称的方式创建模拟。因此对称性不会自发地破坏；当冷热壁之间的温度线性下降时，热对流在这一点上没有出现。在第二阶段，添加随机噪声以触发不稳定性。 这个应用程序在技术上比其他应用程序更先进一点，因为它说明了数据处理器的概念。在本例中，它们用于创建初始条件并触发不稳定性。 此代码中说明的技术： 1. Navier-Stokes和对流-扩散模拟之间的耦合，以使用Boussinesq近似来模拟热流体。 2. 使用数值跟踪器类来决定模拟何时达到稳定状态。 ##### boussinesqThermal3d 限制在热底壁（速度无滑移）和冷顶壁（速度无滑移）之间的流体。侧壁是周期性的。在重力的影响下，形成热对流。热效应通过Boussinesq近似建模：流体是不可压缩的，温度的影响只能通过体力项可见，代表浮力效应。温度场服从对流-扩散方程。 首先以完全对称的方式创建模拟。因此对称性不会自发地破坏；当冷热壁之间的温度线性下降时，热对流在这一点上没有出现。在第二阶段，添加随机噪声以触发不稳定性。 这个应用程序在技术上比其他应用程序更先进一点，因为它说明了数据处理器的概念。在本例中，它们用于创建初始条件并触发不稳定性。 此代码中说明的技术： 1. Navier-Stokes和对流-扩散模拟之间的耦合，以使用Boussinesq近似来模拟热流体。 2. 使用数值跟踪器类来决定模拟何时达到稳定状态。 ##### breakingDam3d 在此示例中，实现了模拟溃坝问题的自由表面流，其中水与障碍物的下游碰撞。该示例说明： 1. 如何建立自由表面流动问题（忽略一相影响的两相问题）。 2. 如何描述自由表面情况下的流动几何。 3. 如何导出数据，特别是如何编写描述自由表面形状的STL文件。 ##### cavity2d 顶盖驱动的二维空腔中流动。空腔是方形的并且没有滑动壁，除了顶壁以恒定速度向右驱动。由于角节点上的速度不连续性，该算例具有挑战性。另一方面，代码非常简单。例如，它可以用作第一个示例，以熟悉Palabos。 此代码中说明的技术： 1. 使用计时器对象对代码进行基准测试。 2. 数据分析：速度、速度范数和涡度的计算。 3. 从 2D 模拟数据创建 GIF 图像和 VTK 文件。 ##### cavity3d 顶盖驱动的三维空腔中流动。在这个三维腔体模拟中，顶盖在平行于两条对角线之一的方向上以恒定速度驱动。由于角节点上的速度不连续性，该算例具有挑战性。另一方面，代码非常简单。例如，它可以用作第一个三维示例，以熟悉Palabos。 此代码中说明的技术： 1. 使用计时器对象对代码进行基准测试。 2. 数据分析：速度、速度范数和涡度的计算。 ##### collidingBubbles3d 以三维形式相互投射的两个气泡之间的碰撞。该应用程序使用`He/Lee`多相流模型。 此代码中说明的技术： 1. 晶格与各种标量场和张量场之间的耦合，使用多相流体的`He/Lee`算法。 2. 通过数据处理函数直接初始化标量场和张量场值。 ##### cylinder2d 围绕通道内的二维圆柱体流动，创建卡门涡街。此示例使用反弹节点来描述圆柱体的形状。出口通过Neumann（零速度梯度）条件建模。 此代码中说明的技术： 1. 使用域函数根据分析处方在特定位置实例化反弹节点。 2. 出口的Neumann速度条件的实例化。 3. 无需计算即可访问内部统计数据以计算平均动能和平均密度等量。 ##### multiComponent2d Rayleigh-Taylor不稳定性发生在二维问题重流体最初位于轻流体之上。该应用程序利用了多相流体的`Shan/Chen`多组分模型。 此代码中说明的技术： 1. 使用Shan/Chen算法在两个Navier-Stokes格子之间进行耦合，以模拟不混溶的多组分流体。 2. 使用单格索引泛函轻松创建空间相关边界条件。 ##### multiComponent3d Rayleigh-Taylor不稳定性发生在重流体最初位于轻流体顶部的三维问题中。该应用程序利用了多相流体的`Shan/Chen`多组分模型。 1. 使用Shan/Chen算法在两个Navier-Stokes格子之间进行耦合，以模拟不混溶的多组分流体。 2. 使用单格索引泛函轻松创建空间相关边界条件。 ##### poiseuille 实现固定的、压力驱动的二维通道流，并与分析Poiseuille曲线进行比较。速度被初始化为零，并且仅缓慢地收敛到预期的抛物线。此应用程序说明了CFD应用程序中的完整生产周期，范围从创建几何图形和定义程序执行过程中的边界条件到评估结果和生成瞬时图形快照。从技术的角度来看，这个展示并不是微不足道的：它实现了例如混合速度/压力边界，并使用分析剖面来设置边界和初始条件，并计算误差。作为第一个Palabos示例，您可能更喜欢查看更简单的代码，例如cavity2d。 此代码中说明的技术： 1. 定义函数（或函数对象）以创建空间相关的初始条件和边界条件。 2. 计算的速度场与解析公式的比较，并计算均方根误差。 3. 定义速度或压力边界，或混合情况，在边界的某些部分使用速度条件，而在其他部分使用压力条件。 ##### rectangularChannel3d 在具有矩形横截面的通道中将泊肃叶流扩展到三维。该流动是静止的，并将数值结果与解析公式进行比较。 ##### womersley 脉动、力驱动的二维通道流的实现，并与分析`Womersley`轮廓进行比较。外力随时间变化为一个简单的余弦；由此产生的速度分布远非微不足道，并且由`Womersley`解决方案描述。 此代码中说明的技术： 1. 外力场的使用（可能与空间有关；在本例中，它与空间无关，但与时间有关）。 2. 在一个方向上实施周期性边界。 #### 目录`examples/codesByTopic` ##### particlesInTube 管道中的三维定常流模拟。被动标量粒子被注入，并写入ASCII文件以进行可视化。 ##### bounceBack 此目录中说明的技术： 1. computeDrag.cpp：通过评估动量交换计算作用在反弹障碍物上的阻力和升力。 ##### boundaryCondition 此目录中说明的技术： 1. neumannOutlets.cpp：实现Neumann型出口的两种不同方式。 * 通过复制前一个站点的速度来获得零速度梯度，并将其用于 Dirichlet 边界条件。 * 复制前一个站点的所有未知粒子群。 ##### dataAnalysis 此目录中说明的技术： 1. cavity3d.cpp：将（有效的）内部统计数据与手动计算的值进行比较，并得出它们相等的结论。 ##### dotList 1. cylinder2d.cpp：使用基于点列表的数据处理器而不是通常的盒装数据处理器来创建复杂域。 ##### couplings 此目录中说明的技术： 1. CoupleVelocityField.cpp：编写一个数据处理器来创建一个标量场和一个块点阵之间的耦合。在此示例中，标量场的值用于设置块格中的边界条件。 ##### include 该目录汇总了各种示例代码使用的常用功能。 ##### io 此目录中说明的技术： 1. checkPointing.cpp：保存模拟状态，以便在程序中断后继续进行。 2. loadGeometry.cpp：从ASCII文件中读取布尔掩码，并根据该布尔掩码的值设置反弹节点。这种方法可用于复杂的几何形状，例如多孔介质。 3. useIOstream.cpp：使用输出流将模拟结果以格式化的方式写入文件或终端。 4. userInput.cpp：从命令行或输入文件获取输入参数。 ##### multiBlock 此目录中说明的技术： 1. manualBlockCavity2d.cpp：使用多块晶格构造器的所有参数来显式地创建块分布。 2. manualBlockCavity3d.cpp：使用多块晶格构造器的所有参数来显式地创建块分布。 3. manualBlockPoiseuille2d.cpp：使用多块点阵构造函数的所有参数来显式地创建块分布。 ##### navierStokesModels 此目录中说明的技术： 1. allModels2d.cpp：列出所有可用于不可压缩Navier-Stokes方程的模型：BGK、正则化、MRT、熵。 ##### nonNewtonian 此目录中说明的技术： 1. carreauPoiseuille.cpp：在二维通道中使用非牛顿Carreau模型。 ##### scalarField 此目录中说明的技术： 1. scalarField2d.cpp：对二维标量域执行基于数组的操作。在此示例中，标量场的值被初始化为与空间相关的正弦波。 2. scalarField3d.cpp：对三维标量场执行基于数组的操作。在这个例子中，标量场的值是一个空间相关的正弦波。 ##### smagorinskyModel 此目录中说明的技术： 1. smagorinskyCavity3D.cpp：具有静态Smaogrinsky模型的顶盖驱动型腔中的LES仿真。请注意，结果不可信，因为边界层没有正确建模。 ### 附录B 部分功能/类的参考 #### 用于模拟设置和其他目的的可变（就地）操作 本附录介绍了块格和数据字段的预定义可变操作。可变意味着原始块格或数据字段被修改。对于数据的分析和提取，通常更适合使用附录用于数据分析和其他目的的非可变操作中介绍的方法。 其中一些函数在其参数列表中接受所谓的“函数”，或行为类似于函数的用户定义对象。这些函数继承自以下三个虚拟基类之一： **OneCellFunctional3D**:此函数用于指定要在格子的子集上执行的完全本地操作。该操作是通过覆盖虚函数`void OneCellFunctional3D::execute(Cell& cell) const`来定义的。 **OneCellIndexedFunctional3D**:有了这个函数，一个完全局部的，但空间相关的动作是在一个格子的一个子集上执行的。动作是通过重写虚函数`void OneCellIndexedFunctional3D::execute(plint iX, plint iY, plint iZ, Cell& cell) const`实现的。 **DomainFunctional3D**:此函数用于定义格子上执行动作的区域。区域是通过覆盖虚函数`bool DomainFunctional3D::operator()(plint iX, plint iY, plint iZ) const`来定义的。 最后，一些函数使用静态函数，由于模板机制，这些函数不受继承层次结构的限制。它们的语法取决于上下文。例如，函数`setBoundaryDensity`的函数`Function f`必须定义一个`T operator() (plint iX, plint iY, plint iZ)`形式的方法，该方法返回每个空间点的密度值。 所有操作都应用于块的子域，当通过`DomainXD`类型的参数`domain`指定时，它是矩形的，或者当通过`DomainFunctional3D`类型的参数`domainFunctional`指定时，它是一般形状的。 后跟星号(*)的参数按指针传递，而所有其他参数均按值或按引用传递。 ##### 块的操作 **apply(lattice, domain, oneCellFunctional\*)**：对域的所有单元应用相同的操作。 **applyIndexed(lattice, domain, oneCellIndexedFunctional\*)**：将依赖于单元格的操作应用于域的所有单元格。 **defineDynamics(lattice, domain, dynamics\*)**：为域的所有单元分配一个新的动力学对象。每个单元格都获得动力学对象的独立副本。 **defineDynamics(lattice, boundingBox, domainFunctional\*, dynamics\*)**：为域的所有单元分配一个新的动力学对象。每个单元格都获得动力学对象的独立副本。参数`boundingBox`是`BoxXD`类型，它必须包含由`domainFunctional`描述的域。出于效率原因，以避免不必要地经常评估域功能。 **defineDynamics(lattice, dotList, dynamics\*)**：为点列表中包含的所有单元分配一个新的动力学对象。每个单元格都获得动力学对象的独立副本。 **defineDynamics(lattice, boolMask, domain, dynamics\*, flag)**:将新的动态对象分配给`boolMask`评估为标记的所有单元格。每个单元格都获得动态对象的独立副本。标志`flag`是布尔值，`boolMask`是`MultiScalaFieldXD`类型，但被评估为布尔值（小于0.5表示false）。如果未指定参数域，则使用整个晶格域。 **defineDynamics(lattice, iX, iY, iZ, dynamics\*)**:将新的动力学对象分配给位置(iX,iY,iZ)的单元格。 **setBoundaryVelocity(lattice, domain, velocity)**:将相同的速度分配给域内为速度实现`Dirichlet`边界条件的所有单元。在所有其他单元格上，此功能无效。 **setBoundaryVelocity(lattice, domain, Function f)**:将与单元相关的速度分配给域内的所有单元，这些单元实现了速度的Dirichlet边界条件。在所有其他单元格上，此功能无效。 **setBoundaryDensity(lattice, domain, rho)**:将依赖于单元的密度分配给域内的所有单元，这些单元实现了密度（或压力）的Dirichlet边界条件。在所有其他单元格上，此功能无效。 **setBoundaryDensity(lattice, domain, Function f)**:将依赖于单元的密度分配给域内的所有单元，这些单元实现了密度（或压力）的Dirichlet边界条件。在所有其他单元格上，此功能无效。 **initializeAtEquilibrium(lattice, domain, rho, velocity)**:以平衡分布初始化域的单元，对所有单元使用相同的密度和速度。平衡态的确切形式由此时驻留在单元上的动力学对象定义。 **initializeAtEquilibrium(lattice, domain, Function f)**:使用与单元相关的密度和速度，以平衡态分布初始化域的单元。平衡态的确切形式由此时驻留在单元上的动力学对象定义。 **stripeOffDensityOffset(lattice, domain, deltaRho)**:将粒子分布函数分解为宏观变量和非平衡部分。然后，从密度中减去值`deltaRho`，并重组分布函数。 **setCompositeDynamics(lattice, domain, compositeDynamics\*)**:使用当前驻留在单元上的动力学作为基础动力学，将复合动力学对象分配给域的所有单元。 **setExternalScalar(lattice, domain, whichScalar, externalScalar)**:为域的所有单元格上的外部标量`whichScalar`分配相同的常数值。 **setExternalVector(lattice, domain, vectorStartsAt, externalVector)**:这里，`externalVector`的类型是`Array`。分配此数组的内容以将值分配给域的每个单元格上的2(2D)或3(3D)外部标量。该函数通常用于在强制格子Boltzmann模拟中初始化外力项。 ##### 标量场的操作 **setToConstant(field, domain, value)**:将域中的所有标量初始化为一个常数值。 **setToFunction(field, domain, Function f)**:将域中的所有标量初始化为一个空间函数值（函数f采用两个(2D)或三个(3D)整数坐标值，并返回T）。 **setToCoordinate(field, domain, index)**:将域中的每个标量替换为其空间位置的索引分量（x=0、y=1 和 z=2 表示索引）。 **apply(Function f, field, domain)**:将任何函数或函数f（接受T并返回T）应用于域。 ##### 标量场操作：A = A op alpha **addInPlace(field, double scalar, Box2D domain)**:将标量值标量添加到域内字段的每个值。 **subtractInPlace(field, double scalar, Box2D domain)**:从域内字段的每个值中减去标量值标量。 **multiplyInPlace(field, double scalar, Box2D domain)**:将字段的每个值乘以域内的标量标量。 **divideInPlace(field, double scalar, Box2D domain)**:将字段的每个值除以域内的标量标量。 ##### 标量场操作:A = A op B **addInPlace(A, B, domain)**:在domain中计算$A = A + B$ **subtractInPlace(A, B, domain)**:在domain中计算$A = A - B$ **multiplyInPlace(A, B, domain)**:在domain中计算$A = A \times B$ **divideInPlace(A, B, domain)**:在domain中计算$A = A / B$ ##### 张量场操作 **setToConstant(field, domain, value)**:将域中的所有张量/向量初始化为一个常数值（值的类型为 `Array`）。 **setToFunction(field, domain, Function f)**:将域中的所有张量/向量初始化为一个空间函数值（函数f采用两个(2D)或三个(3D)整数坐标值，并返回`Array`）。 **setToCoordinates(field, domain, index)**:将域中的每个向量替换为其空间位置。张量场必须包含维度为2(2D)或3(3D)的向量。 ##### 张量场操作:A = A op B **addInPlace(A, B, domain)**:在domain中计算$A = A + B$ **subtractInPlace(A, B, domain)**:在domain中计算$A = A - B$ **multiplyInPlace(A, B, domain)**:在domain中计算$A = A \times B$ **divideInPlace(A, B, domain)**:在domain中计算$A = A / B$ #### 用于数据分析和其他目的的非可变操作 本附录介绍了用于块格和数据字段的预定义不可变操作。不可变意味着原始块格或数据字段不被修改，并且结果存储在新的块格或数据字段中。这些操作通常用于数据后处理，因此以“数据分析”的名义呈现。对于模拟期间的实际计算，使用附录用于模拟设置和其他目的的可变（就地）操作小节中的介绍通常更合适。 运算符在内部被定义为数据处理器，并通过便利函数呈现在用户界面中，如下所示的函数`computeDensity()`： **std::unique_ptr > computeDensity(MultiBlockLattice3D& lattice)**:此函数将请求的操作应用于整个格域，并构造一个新的标量域来存储结果。每个操作者都可以使用这种形式的功能。 **std::unique_ptr > computeDensity(MultiBlockLattice3D& lattice, domain)**:在这种情况下，该操作仅适用于子域，并且生成的标量域具有与域对应的大小。这种形式的函数适用于每个操作符，并且之后会用该形式来说明运算符的使用。 **void computeDensity(MultiBlockLattice3D& lattice, MultiScalarField3D& density,domain)**:在这个函数中，没有分配新的字段，因为结果存储在字段密度中，由参数给这个函数。这种形式的函数始终可用，但返回标量值的运算符除外，例如`computeAverageDensity()`。 更多关于使用非可变操作的解释可以在章节数据评估中找到，函数的确切语法在`Palabos`参考指南中给出。 ##### 块->标量 **computeAverageDensity(lattice, domain)**：返回在域上计算的密度平均值。 **computeAverageRhoBar(lattice, domain)**：返回在域上计算的rhoBar的平均值。请注意，虽然从概念上讲，命令`computeAverageDensity()`等效于`Descriptor::fullRho(computeRhoBar())`，但后者可能较少受到舍入误差的影响。 **computeAverageEnergy(lattice, domain)**：返回速度的范数平方的平均值，并除以2。 ##### 块->块 **extractSubDomain(lattice, domain)**：从格子的域中复制值并将它们存放在相应大小的新格子中。 ##### 块->标量场/张量场 **computeDensity(lattice, domain)**：计算每一个单元的$\rho$。 **computeRhoBar(lattice, domain)**：计算每一个单元的$\bar{\rho}$。 **computeKineticEnergy(lattice, domain)**：计算每一个单元的动能。 **computeVelocityNorm(lattice, domain)**：计算每一个单元的速度范数。 **computeVelocityComponent(lattice, domain, iComponent)**：计算每一个单元的速度分量（对于二维问题来说，iComponent的取值范围为0和1；对于三维问题来说，iComponent的取值范围为0，1，2）。 **computeVelocity(lattice, domain)**：计算每个单元上速度的所有分量。 **computeDeviatoricStress(lattice, domain)**：计算每个单元上应力张量$\Pi$的非平衡部分。结果取决于模型，由在单元位置定义的动力学对象确定。结果存储在3分量(2D)或6分量(3D)张量场中（分量数量减少，因为张量是对称的）。要知道特定分量的索引，请使用`SymmetricTensorImpl::xz`之类的符号来表示三维对称张量的xy分量。 **computeStrainRateFromStress(lattice, domain)**：计算应变率$S=1/2(\nabla{U}+\nabla{U}^T)$，假设它与偏应力成正比。结果取决于模型，并由在单元位置定义的动力学对象确定。结果存储在3分量(2D)或6分量(3D)张量场中（分量数量减少，因为张量是对称的）。要知道特定分量的索引，请使用`SymmetricTensorImpl::xz`之类的符号表示三维对称张量的xy分量。 **computePopulation(lattice, domain, iPop)**：提取每个单元格上的分布函数数量iPop，并将结果存放在标量场中。iPop的值范围从0到q-1。 ##### 标量场->标量 **computeAverage(scalarField, domain)**：返回域中标量字段中值的平均值。 **computeBoundedAverage(scalarField, domain)**：使用梯形规则以二阶精度近似域上的积分（与`computeAverage()`相比，边界节点的权重小于块节点。 **computeMin(scalarField, domain)**：在域上返回标量字段中的最小值（不返回此值的位置）。 **computeMax(scalarField, domain)**：在域上返回标量字段中的最大值（不返回此值的位置）。 ##### 标量场->标量场 **extractSubDomain(scalarField, domain)**：从`scalarField`的域中复制值并将它们存放在相应大小的新标量域中。 ##### 标量场：B = A op alpha **add(scalar, field, domain)**:返回在域的每个元素上计算`scalar + field(x,y,z)`的结果。 **add(field, scalar, domain)**:返回在域的每个元素上计算`field(x,y,z) + scalar`的结果。 **subtract(scalar, field, domain)**:返回在域的每个元素上计算`scalar - field(x,y,z)`的结果。 **subtract(field, scalar, domain)**:返回在域的每个元素上计算`field(x,y,z) - scalar`的结果。 **multiply(scalar, field, domain)**:返回在域的每个元素上计算`scalar * field(x,y,z)`的结果。 **multiply(field, scalar, domain)**:返回在域的每个元素上计算`field(x,y,z) * scalar`的结果。 **divide(scalar, field, domain)**:返回在域的每个元素上计算`scalar / field(x,y,z)`的结果。 **divide(field, scalar, domain)**:返回在域的每个元素上计算`field(x,y,z) / scalar`的结果。 ##### 标量场->标量场：C = A op B **add(A, B, domain)**:返回C = A + B，域中的每个单元都需要计算。 **subtract(A, B, domain)**:返回C = A - B，域中的每个单元都需要计算。 **multiply(A, B, domain)**:返回C = A * B，域中的每个单元都需要计算。 **divide(A, B, domain)**:返回C = A / B，域中的每个单元都需要计算。 ##### 张量场->张量场 **extractSubDomain(tensorField, domain)**:从`tensorField`的域中复制值并将它们存放在相应大小的新张量域中。 **extractComponent(tensorField, domain, iComponent)**:将组件`iComponent`从`tensorField`的每个单元格复制到生成的标量场中。 **computeVorticity(velocity, domain) [3D only]**:计算三维速度场中每个单元的涡量，并将结果存储在三维涡量场中。涡量是通过在大部分域上的中心（二阶精确）有限差分格式以及通过边界节点上的不对称最近邻格式来评估的。 **computeBulkVorticity(velocity, domain) [3D only]**:计算三维速度场中每个单元的涡量，并将结果存储在三维涡度场中。涡量通过对所有单元的中心（二阶精确）有限差分方案进行评估。如果域包括模拟的非周期性边界节点，您应该使用`computeVorticity()`来考虑边界效应。 **computeStrainRate(velocity, domain)**:计算速度场每个单元的应变率$S=1/2(\nabla{U}+\nabla{U}^T)$，并将结果存储在3分量(2D)或6分量(3D)中 张量场（由于张量是对称的，因此减少了组件的数量。要知道特定组件的索引，请使用`SymmetricTensorImpl::xz`之类的符号表示三维对称张量的xy组件。应变率是通过在大部分域上的中心（二阶精确）有限差分方案以及通过边界节点上的非对称最近邻方案来评估的。 **computeBulkStrainRate(velocity, domain)**:计算速度场每个单元的应变率，并将结果存储在3分量(2D)或6分量(3D)张量场中。应变率通过对所有单元的中心（二阶精确）有限差分方案进行评估。如果域包括模拟的非周期性边界节点，您应该使用`computeStrainRate()`来考虑边界效应。 ##### 张量场->标量场 **computeNorm(tensorField, domain)**:计算`tensorField`中每个单元的欧几里得范数。 **computeNormSqr(tensorField, domain)**:计算`tensorField`中每个单元的欧几里得范数的平方。 **computeSymmetricTensorNorm(tensorField, domain)**:计算对称张量`tensorField`中每个单元的欧几里得范数。由于对称性，假设`tensorField`仅存储一半的非对角线分量。 **computeSymmetricTensorNormSqr(tensorField, domain)**:计算对称张量`tensorField`中每个单元的欧几里得范数的平方。由于对称性，假设`tensorField`仅存储一半的非对角线分量。 **computeSymmetricTensorTrace(tensorField, domain)**:计算对称张量`tensorField`中每个单元的轨迹。由于对称性，假设`tensorField`仅存储一半的非对角线分量。 **computeVorticity(velocity, domain) [2D only]**:计算二维速度场中每个单元的涡量，并将结果存储在标量场中。涡量是通过在大部分域上的中心（二阶精确）有限差分格式以及通过边界节点上的不对称最近邻格式来评估的。 **computeBulkVorticity(velocity, domain) [2D only]**:计算二维速度场中每个单元的涡量，并将结果存储在标量场中。涡量通过对所有单元的中心（二阶精确）有限差分方案进行评估。如果域包括模拟的非周期性边界节点，您应该使用`computeVorticity()`来考虑边界效应。 ##### 张量场：C = A op B **add(A, B, domain)**: 返回C = A + B，计算域中的每个单元。 **subtract(A, B, domain)**: 返回C = A - B，计算域中的每个单元。 **multiply(A, B, domain)**: 返回C = A * B，计算域中的每个单元。 **divide(A, B, domain)**: 返回C = A / B，计算域中的每个单元。 ### 附录C 版权和许可协议 ## 参考文献