LBM/PALABOS：PALABOS介绍与安装

[toc]

软件简介

Palabos 库是通用计算流体动力学 (CFD) 的框架，其内核基于格子玻尔兹曼 (LB) 方法。它既可用作研究工具，也可用作工程工具：其编程界面简单明了，可以相对轻松地设置流体流动模拟，或者，如果您了解格子玻尔兹曼方法，则可以使用自己的方法扩展该库楷模。 Palabos 代表并行格子玻尔兹曼求解器。
在线手册、在线教程

什么是LBM？

概况

晶格玻尔兹曼方法是计算流体力学中的一种现代方法。它通常用于求解不可压缩的、时变的Navier-Stokes方程。
它的优势在于能够很容易地表示复杂的物理现象，从多相流动到流体与周围环境之间的化学相互作用。
该方法起源于对流体的分子描述，并且可以直接纳入来自分子间相互作用知识的物理术语。
因此，它在基础研究中是一个宝贵的工具，因为它使理论的阐述和相应数值模型的制定之间的周期缩短。同时，它已被证明是一种高效和方便的替代传统的解决方案，为各种各样的工业问题。

起源

晶格玻尔兹曼通常被简单地认为是玻尔兹曼方程的数值解算器。玻尔兹曼方程是分子水平上的纳维-斯托克斯方程的类似物，它描述了一个被称为概率分布函数的统计量的时空动力学，它被定义在6维相空间中。
在分子水平的描述上，模型所涵盖的物理现象的数量比在Navier-Stokes方程的水动力水平上要大。这是因为玻尔兹曼方程不受时间尺度分离的影响，能够描述具有大分子平均自由路径的非水动力状态下的流体。
此外，该分子模型能够捕捉到摩擦、扩散和温度输运等输运现象，并推导出相应的输运系数。
从历史上看，晶格玻尔兹曼模型是从所谓的元胞自动机模型演变而来的，独立于上面展示的理论。元胞自动机（Cellular Automaton）是用整数值描述的离散状态演化的计算机模型(与晶格玻尔兹曼变量的浮点表示相反)。
在我们的例子中，这种状态代表了相互作用的伪分子，流体的“介观”粒子的位置和速度。为了更好地理解这种模型的连续介质物理学，可以将元胞自动机的玻尔兹曼方程写下来。这个方程作用于实值量，但它描述了离散相空间中的一些动力学，在下面我们称之为晶格。晶格玻尔兹曼模型是在首次尝试短路元胞自动机过程并直接模拟相应离散玻尔兹曼方程的动力学时诞生的。
虽然后来人们了解到晶格玻尔兹曼方程可以直接从连续玻尔兹曼方程推导出来，但记住它的真正起源是很有趣的。事实上，这揭示了其特别优雅的公式的一些原因，使相应的计算机程序能够直接实现。

为什么使用晶格玻尔兹曼?

每个人都同意晶格玻尔兹曼在理论、代码实现和应用程序的外观和感觉层面上与经典CFD不同。但是为什么你会选择晶格玻尔兹曼而不是其他方法呢?
在这一点上意见分歧更大一些。然而，从我们发展Palabos的经验来看，我们认为晶格玻尔兹曼的根本优势是效率。该方法从一开始就设计为在高性能硬件上运行，并且它适应复杂的物理或复杂的算法。
当效率允许解决以前无法解决的问题，或者只是不够精确的问题时，就会带来全新的用户体验。
像PALABOS这样的晶格玻尔兹曼CODE的典型成就是
- 数据预处理和网格生成的时间只占整个模拟的一小部分。
- 并行数据分析，后处理和评估。
- 完全解决多相流与小液滴和气泡。
- 完全解决流动通过复杂的几何形状和多孔介质。
- 具有传热和化学反应的复杂、耦合流动。
晶格玻尔兹曼在建模方法的层面上也提供了优势。该方法在某种意义上是混合的，因为求解器是基于网格的(如CFD中通常的那样)，但也继承了基于粒子的方法的某些方面，因为它起源于细胞自动机。因此，很容易将其与嵌入粒子方法相结合，用于沉积物输运或其他现象。

晶格玻尔兹曼和高性能计算

与经典CFD相比，晶格玻尔兹曼乍一看似乎相当消耗资源:该模型中使用的离散概率分布函数比Navier-Stokes方程的流体动力变量(每个节点9个实值，而二维不可压缩解算器为3个)需要更多的存储空间。
然而，这从来都不是一个真正的问题，特别是在现代计算机上，它被出色的计算效率大大补偿了。由于其明确的公式和精确的平流算子，晶格玻尔兹曼格式只涉及非常有限数量的每个计算节点的浮点运算。
此外，晶格玻尔兹曼方法特别适合于并行体系结构上的计算，即使具有缓慢的互连网络。这种优势也扩展到其他类型的高性能硬件，如通用图形处理单元(gpgpu)，这是晶格玻尔兹曼实现的目标。

下载安装

LINUX（windows子系统Ubuntu20.04测试）安装：

# 安装必要依赖
sudo apt install gcc clang clang-format cmake make
# 安装可选依赖，已有mpi 这里不安装
sudo apt install imagemagick libhdf5-dev libhdf5-mpi-dev ccache
# 从gitlab clone palabos包
git clone https://gitlab.com/unigespc/palabos.git
# 或者下载、解压
wget https://gitlab.com/unigespc/palabos/-/archive/v2.3.0/palabos-v2.3.0.tar.gz
tar xzvf palabos-v2.3.0.tar.gz

测试

测试cavity2d例子

cd palabos/examples/showCases/cavity2d/build/
cmake ..
make
cd ..
./cavity2d
# 并行计算
mpirun -np 4 ./cavity2d

uNorm128000

预处理和后处理软件（开源）

ParaView ：Paraview为VTK库提供了一个图形用户界面，是一个很好的科学数据可视化工具。Palabos模拟的结果可以保存在一个VTK格式，然后用Paraview处理。
Octave：这个命令行解释器与Matlab非常相似。它对于处理由Palabos生成的ASCII数据非常有用，例如创建图形。
1
sudo apt-get install octave
ImageMagick：不需要外部库，Palabos就可以生成PPM格式的快照图像。如果安装了ImageMagick，这些图片可以自动转换成更常见的GIF格式。此外，ImageMagick可用于从一系列GIF图像中生成动画GIF。
Meshlab：Meshlab将表面网格文件从许多不同的格式转换为STL格式，Palabos需要它来设置几何图形。此外，Meshlab通常能够纠正STL文件中的错误，例如，删除零面积三角形或恢复表面法线。

教程

Tutorial 1

为了说明Palabos中的编程概念，模拟了一个周期性的边界的没有边界条件的初值问题的时间演化；初始条件为零速度和恒定密度；在场的平方子域上，密度初始化为稍高的值，以产生扰动，然后产生流型。

应当特别指出，初始条件并不代表不可压缩流体的状态，因为零速度与非零密度是不相容的。

经过足够的迭代，流动自动收敛到与不可压缩或略可压缩流的物理相容的情况。

cd ~/palabos/examples/tutorial/tutorial_1
cd build
cmake ..
make
./tutorial_1_1
# Writing GIF file at iT=0
# Writing GIF file at iT=40
# Writing GIF file at iT=80
# Writing GIF file at iT=120
# Writing GIF file at iT=160
# Writing GIF file at iT=200
# ...

merged

代码解读


/* Code 1.1 in the Palabos tutorial
 */

#include <iomanip>
#include <iostream>

#include "palabos2D.h" // 允许访问Palabos项目中的所有声明
#include "palabos2D.hh" // 为了保证访问完整的模板代码，例如，使用不同于标准数据类型double的数据类型实例化模拟

using namespace plb; // 保证了对Palabos代码中的名称的自动访问
using namespace std;

typedef double T; // 将在二维D2Q9晶格上使用双精度浮点数执行。
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor

// Initialize the lattice at zero velocity and constant density, except
//   for a slight density excess on a square sub-domain.
void defineInitialDensityAtCenter(MultiBlockLattice2D<T, DESCRIPTOR> &lattice)
{
    // The lattice is of size nx-by-ny
    const plint nx = lattice.getNx();
    const plint ny = lattice.getNy();

    // Create a Box2D which describes the location of cells with a slightly
    //   higher density. 子域高密度初始化
    plint centralSquareRadius = nx / 6;
    plint centerX = nx / 3;
    plint centerY = ny / 4;
    Box2D centralSquare(
        centerX - centralSquareRadius, centerX + centralSquareRadius, centerY - centralSquareRadius,
        centerY + centralSquareRadius);

    // All cells have initially density rho ...初始化所有晶格单元在密度恒定、速度为零的平衡分布
    T rho0 = 1.;
    // .. except for those in the box "centralSquare" which have density
    //    rho+deltaRho
    T deltaRho = 1.e-4;
    Array<T, 2> u0((T)0, (T)0);

    // Initialize constant density everywhere.
    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), rho0, u0);
	// 重新初始化子域高密度
    // And slightly higher density in the central box.
    initializeAtEquilibrium(lattice, centralSquare, rho0 + deltaRho, u0);

    lattice.initialize(); // 初始化操作后调用initialize方法，为模拟准备晶格。
}
// 模拟过程主函数
int main(int argc, char *argv[])
{
    plbInit(&argc, &argv); // 函数plbInit必须在程序的最开始强制调用，以保证顺序和并行程序运行之间的结果一致。
    global::directories().setOutputDir("./tmp/"); // 指定输出文件夹

    const plint maxIter = 1000;  // Iterate during 1000 steps.
    const plint nx = 600;        // Choice of lattice dimensions.
    const plint ny = 600;
    const T omega = 1.;  // Choice of the relaxation parameter
	// 在创建多块期间，默认动力学(在本例中为BGK)	
    MultiBlockLattice2D<T, DESCRIPTOR> lattice(nx, ny, new BGKdynamics<T, DESCRIPTOR>(omega));
	// 在Palabos中，周期边界以一种特殊的方式处理。它们只能通过调用方法periodicperiodic()强加于块晶格的外部边界。
    lattice.periodicity().toggleAll(true);  // Use periodic boundaries.

    defineInitialDensityAtCenter(lattice);

    // Main loop over time iterations.
    for (plint iT = 0; iT < maxIter; ++iT) {
        if (iT % 40 == 0) {  // Write an image every 40th time step.
            pcout << "Writing GIF file at iT=" << iT << endl;
            // Instantiate an image writer with the color map "leeloo".
            ImageWriter<T> imageWriter("leeloo");
            // Write a GIF file with colors rescaled to the range of values
            //   in the matrix
            imageWriter.writeScaledGif(createFileName("u", iT, 6), *computeVelocityNorm(lattice));
        }
        // Execute lattice Boltzmann iteration.
        lattice.collideAndStream();
    }
}

Tutorial 2

修改CMakeList.txt内容：
1
set(EXECUTABLE_NAME "tutorial_1_2")
编译与运行：
1
2
3
4
cd build
cmake ..
make
./tutorial_1_2
结果：

与第一个例子的区别：

const plint maxIter = 1000;  // Iterate during 1000 steps.
const plint nx = 600;        // Choice of lattice dimensions.
const plint ny = 600;
const T omega = 1.;  // Choice of the relaxation parameter

T rho0 = 1.;  // All cells have initially density rho ...
// .. except for those inside the disk which have density
//    rho+deltaRho
T deltaRho = 1.e-4;
Array<T, 2> u0(0, 0);
// 为子域的每个单元分配恒定的密度，与之前的initializeAtEquilibrium具有相同效果
void initializeConstRho(plint, plint, T &rho, Array<T, 2> &u)
{
    u = u0;
    rho = rho0 + deltaRho;
}
// 高密度区域改成圆盘
void initializeRhoOnDisk(plint iX, plint iY, T &rho, Array<T, 2> &u)
{
    plint radius = nx / 6;
    plint centerX = nx / 3;
    plint centerY = ny / 4;
    u = u0;
    if ((iX - centerX) * (iX - centerX) + (iY - centerY) * (iY - centerY) < radius * radius) {
        rho = rho0 + deltaRho;
    } else {
        rho = rho0;
    }
}

// Initialize the lattice at zero velocity and constant density, except
//   for a slight density excess on a circular sub-domain.
void defineInitialDensityAtCenter(MultiBlockLattice2D<T, DESCRIPTOR> &lattice)
{
    // Initialize constant density everywhere.
    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), rho0, u0);

    // And slightly higher density in the central box.
    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), initializeRhoOnDisk);

    lattice.initialize();
}

Tutorial 3

Compilation options

Tutorial 4

Data analysis

主要在main函数里进行数据分析

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);
    global::directories().setOutputDir("./tmp/");

    MultiBlockLattice2D<T, DESCRIPTOR> lattice(nx, ny, new BGKdynamics<T, DESCRIPTOR>(omega));

    lattice.periodicity().toggleAll(true);  // Set periodic boundaries.

    defineInitialDensityAtCenter(lattice);

    // First part: loop over time iterations.
    for (plint iT = 0; iT < maxIter; ++iT) {
        lattice.collideAndStream();
    }

    // Second part: Data analysis.
    Array<T, 2> velocity;
    lattice.get(nx / 2, ny / 2).computeVelocity(velocity); // 方法get提供对晶格cell的访问。cell是一个对象，通过它不仅可以访问在cell上定义的实际变量，还可以执行有用的计算，例如计算宏观变量。在本例中，速度是在域中间的单元格上计算的
    pcout << "Velocity in the middle of the lattice: (" << velocity[0] << "," << velocity[1] << ")"
          << endl;

    pcout << "Velocity norm along a horizontal line: " << endl;
    Box2D line(0, 100, ny / 2, ny / 2);
    pcout << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl; // 在提取的域上计算速度范数并打印到终端。

    plb_ofstream ofile("profile.dat"); // 写入文件，确保使用数据类型plb_ofstream而不是ofstream来保证并行程序的工作条件
    ofile << setprecision(3) << *computeVelocityNorm(*extractSubDomain(lattice, line)) << endl;

    pcout << "Average density in the domain: " << computeAverageDensity(lattice) << endl; // 计算平均密度
    pcout << "Average energy in the domain: " << computeAverageEnergy(lattice) << endl; // 计算平均动能
}

octave画图

octave # 进入octave
load profile.dat
plot(profile)
xlabel ("Time (s)")
ylabel ("Velocity (m/s)")
print(fig, "velocity.png")

velocity

Tutorial 5

泊肃叶（Poiseuille）流例子：这种流动在无滑移壁面的通道中发展，并且流动速度处处平行于通道壁面

在仿真中，采用解析解实现了区域进出口速度的Dirichlet边界条件。

采用零速度场作为初始条件，开始向域内各点的Poiseuille解收敛。

注意，模拟在初始条件下的零速度和边界条件下的有限速度之间存在不连续，如果在处理参数时模拟在数值上变得不稳定，请尝试增加晶格分辨率或降低晶格单位的速度。

修改CMakeList.txt内容：
1
set(EXECUTABLE_NAME "tutorial_1_5")
编译与运行：
1
2
3
4
cd build
cmake ..
make
./tutorial_1_5
结果：

代码：


#include <iomanip>
#include <iostream>
#include <vector>

#include "palabos2D.h"
#include "palabos2D.hh"

/* Code 1.5 in the Palabos tutorial
 */

using namespace plb;
using namespace std;

typedef double T;
#define DESCRIPTOR plb::descriptors::D2Q9Descriptor

/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const &parameters)
{
    T y = (T)iY / parameters.getResolution();
    return 4. * parameters.getLatticeU() * (y - y * y);
} //这个函数在给定高度的通道中，以给定的最大速度在通道中间计算抛物型泊泽维尔剖面

/// A functional, used to initialize the velocity for the boundary conditions
template <typename T>
class PoiseuilleVelocity {
public:
    PoiseuilleVelocity(IncomprFlowParam<T> parameters_) : parameters(parameters_) { }
    /// This version of the operator returns the velocity only,
    ///    to instantiate the boundary condition.
    void operator()(plint, plint iY, Array<T, 2> &u) const
    {
        u[0] = poiseuilleVelocity(iY, parameters);
        u[1] = T();
    }
    /// This version of the operator returns also a constant value for
    ///    the density, to create the initial condition.
    void operator()(plint, plint iY, T &rho, Array<T, 2> &u) const
    {
        u[0] = poiseuilleVelocity(iY, parameters);
        u[1] = T();
        rho = (T)1;
    }

private:
    IncomprFlowParam<T> parameters;
};

void channelSetup(
    MultiBlockLattice2D<T, DESCRIPTOR> &lattice, IncomprFlowParam<T> const &parameters,
    OnLatticeBoundaryCondition2D<T, DESCRIPTOR> &boundaryCondition)
{
    // Create Velocity boundary conditions.
    boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);

    // Specify the boundary velocity.
    setBoundaryVelocity(lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters));

    // Create the initial condition.
    initializeAtEquilibrium(lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters));

    lattice.initialize();
}

void writeGifs(MultiBlockLattice2D<T, DESCRIPTOR> &lattice, plint iter)
{
    const plint imSize = 600;
    ImageWriter<T> imageWriter("leeloo");
    imageWriter.writeScaledGif(
        createFileName("u", iter, 6), *computeVelocityNorm(lattice), imSize, imSize);
}

int main(int argc, char *argv[])
{
    plbInit(&argc, &argv);

    global::directories().setOutputDir("./tmp/");

    // Use the class IncomprFlowParam to convert from
    //   dimensionless variables to lattice units, in the
    //   context of incompressible flows.
    IncomprFlowParam<T> parameters(
        (T)1e-2,  // Reference velocity (the maximum velocity
                  //   in the Poiseuille profile) in lattice units.
        (T)100.,  // Reynolds number
        100,      // Resolution of the reference length (channel height).
        2.,       // Channel length in dimensionless variables
        1.        // Channel height in dimensionless variables
    );
    const T imSave = (T)0.1;  // Time intervals at which to save GIF
                              //   images, in dimensionless time units.
    const T maxT = (T)3.1;    // Total simulation time, in dimensionless
                              //   time units.

    writeLogFile(parameters, "Poiseuille flow");

    MultiBlockLattice2D<T, DESCRIPTOR> lattice(
        parameters.getNx(), parameters.getNy(),
        new BGKdynamics<T, DESCRIPTOR>(parameters.getOmega()));

    OnLatticeBoundaryCondition2D<T, DESCRIPTOR> *boundaryCondition =
        createLocalBoundaryCondition2D<T, DESCRIPTOR>();

    channelSetup(lattice, parameters, *boundaryCondition);

    // Main loop over time iterations.
    for (plint iT = 0; iT * parameters.getDeltaT() < maxT; ++iT) {
        if (iT % parameters.nStep(imSave) == 0 && iT > 0) {
            pcout << "Saving Gif at time step " << iT << endl;
            writeGifs(lattice, iT);
        }
        // Execute lattice Boltzmann iteration.
        lattice.collideAndStream();
    }

    delete boundaryCondition;
}

← LAMMPS计算TIP4P/ICE水温度分布时，低于系统平均温度 LBM/PALABOS：PALABOS代码风格与开发扩展方法 →