翻译了一下官方文档的两个部分，从中选取了一些比较重要的记录下来：

1 Palabos程序的时间周期

LBM采用显式求解器的方法：

碰撞：

$$
f_{i}^{\star}(x, t)=f_{i}(x, t)-\frac{\Delta t}{\tau}\left[f_{i}(x, t)-f_{i}^{\mathrm{eq}}(x, t)\right] \quad
$$

迁移：
$$
f_{i}\left(x+c_{i} \Delta t, t+\Delta t\right)=f_{i}^{*}(x, t)
$$

可以看到，在碰撞步没有发生时间推进，得到的是一个中间变量$f_{i}^{\star}(x, t)$。时间推进是在迁移步完成之后进行。

Palabos采用的亦是同样的显式求解方法。

  1. 首先，所有流体变量都在时间t定义。此时粒子群处于预碰撞状态（也称为“流入分布函数”）。
  2. 碰撞操作符适用于所有单元格。它们现在处于碰撞后状态（也称为“流出变量”）。
  3. 流操作符应用于格子。
  4. 执行所有数据处理器，以执行非局部操作或格子之间的耦合。
  5. 分布函数现在再次处于预碰撞状态，但是时间步往前推了一步：变成t+1。

通过调用方法lattice.collide（）执行碰撞步骤（步骤2），并使用方法lattice.stream（）调用迁移步骤（步骤3 ）。这两个方法也可以组合成对lattice.collideAndStream（）的单个调用。在大多数硬件平台上，collideAndStream（）版本在计算上更有效，因为它只通过格子的内存来执行一次。

可以通过调用方法lattice.executeInternalProcessors（）手动执行数据处理器，如执行，集成和包装数据处理功能部分中所述。然而，这种情况很少发生，因为数据处理器也在函数stream（）和函数collideAndStream（）的末尾自动执行。

如果需要在不执行数据处理的情况下调用stream（）或collideAndStream，例如调试程序，可以使用这些函数的域版本，例如lattice.collideAndStream（lattice.getBoundingBox（））。或者是专门指定一个BoxXD：lattice.collideAndStream（new BoxXd(nx,ny)）

在Palabos中，数据处理器总是在碰撞迁移循环之后执行。因此，在分布函数处于预碰撞状态的时刻，即，总是在碰撞迁移循环之后执行非局部操作。

在设置初始条件之后，要进行初始化。这是通过调用方法lattice.initialize（）来实现的。在开始第一个迭代步骤之前。该方法执行一次数据处理器，并在块内执行内部通信步骤，以保证其内部状态一致。

2 评估数据

建议仅在系统处于预碰撞状态（进入群体）时才计算流体动力学变量。虽然这种区别对于保守变量密度和速度并不重要（它们在碰撞前和碰撞后状态相同），但对于非保守变量（如应力张量）非常重要。只有当它们在预碰撞状态下计算时，非守恒速度矩才能与流体动力学变量相关。请注意，如果使用方法collideAndStream（），则不会有出错的风险，因为无论如何都无法访问冲突后变量。

计算平均能量：

1
pcout  <<  computeAverageEnergy （lattice ） <<  endl ;

3 生成图像

可以从整个域，子域或标量场的切片产生图像。默认情况下，Palabos以PPM格式写入图像，这种ASCII格式很容易生成，无需外部图形库。如果安装了包ImageMagick（特别是此包中包含的命令行工具convert），则还可以以更标准的GIF格式写入图像。这在Linux下也可以在Windows下运行。目录examples / showCases中的大多数示例说明了如何编写GIF图像。

  1. 为了写入图像，首先需要定义一个写入图像的类，就是ImageWriter类：
1
ImageWriter imageWriter("leeloo");

参数中的leeoo是一个string类的对象，它的数值并不是随意取的。如下图所示，在内部调用中会根据这个参数的不同而选取不同的绘图方法。它代表用于将标量映射到三RGB颜色方案的颜色映射方案。其数值可以有：earth, water, air, fire, 和leeloo。

map
map
  1. 下一步，则是要生成PPM或者GIF图像了。

可以使用调整到固定颜色值范围的色彩图，或者使用缩放范围来自动调整以适应数据的最小值和最大值。使用GIF格式，图像也可以重新缩放为不同的大小：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// 写入一个PPM图像，数值范围 [minValue, maxValue].
imageWriter.writePpm(scalarField, "myImage", minValue, maxValue);

// 写入一个PPM图像.数值范围自适应选择
imageWriter.writeScaledPpm(scalarField, "myImage");

// 写入一个GIF图像 ，数值范围 [minValue, maxValue].
imageWriter.writeGif(scalarField, "myImage", minValue, maxValue);

// 写入一个GIF图像.数值范围自适应选择.
imageWriter.writeScaledGif(scalarField, "myImage");

// 写入一个PPM图像.数值范围自适应选择，然后重新对其进行缩放
//  到一个 siyeX-by-sizeY 范围内，不过纵横比不变.
imageWriter.writeScaledGif(scalarField, "myImage", sizeX, sizeY);

使用此方法从3D模拟生成快照也非常有用。在这种情况下，需要从3D数据中提取切片，其形式为具有一个拟2维的图（你可以理解为高度只有一层网格）的3D Box：

1
2
Box3D slice(0,nx-1, 0,ny-1, zValue, zValue);
ImageWriter("earth").writeScaledGif(*computeVelocity(lattice, slice));

3 用于后处理的VTK输出

来自块格（lattice）的标量场或张量场的所有数据都可以以VTK数据格式写入文件中。从那里，它可以通过诸如Paraview的科学数据可视化工具进一步后处理。

也就是说，VTK 是一个非常重要的数据输出。因为它可以直接用于后处理。

Palabos中使用的VTK格式是基于Base64格式的数据的无损二进制表示（它是基于ASCII的二进制格式）。使用保持模拟数据的完整数值精度的格式通常是不必要的的，因为后处理操作通常需要比CFD模拟更低的数值精度。

将双精度数据转换为单精度，可以节省下不少内存，如下面的示例所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 评估离散参数, 以便将参数以无量纲的形式写入VTK
T dx = parameters.getDeltaX();
T dt = parameters.getDeltaT();

// 创建一个VtkImageOutput3D的对象，通过dx指定单元格的宽度
// 其中模板参数为T
VtkImageOutput3D vtkOut("simulationData.dat", dx);

// 添加一个3D张量场到VTK文件
// 其值为速度，单位为dx/dt，是一种无量纲的形式
// 强制转换器精度类型为单精度浮点型
vtkOut.writeData<3,float>(*computeVelocity(lattice), "velocity", dx/dt);

// 添加一个3D张量场到VTK文件，表示涡量
vtkOut.writeData<3,float>(*computeVorticity(*computeVelocity(lattice)), "vorticity", 1./dt);

// 添加一个标量场，表示密度
vtkOut.writeData<1,float>(*computeDensity(lattice), "density", 1.);

上述的列子可以看出，你可以在一个VtkImageOutput3D类型的对象中添加多种变量，比如上面的例子添加了速度，涡量还有密度。此外，这些变量的分量也是可以多样的，比如有三个分量的张量场，以及一维的标量场。

当然，对于2D也有类似的接口：VtkImageOutput2D。不过，分析2D的数据，Paraview不一定是最佳的选择，可以直接输出数值，通过MATLAB等工具进行可视化后处理。

5 检查点：保存并加载模拟状态

使用函数saveBinaryBlock可以保存仿真结果。Palabos类型的数据被保存下来，可以重启动仿真分析。

1
saveBinaryBlock(lattice, "checkpoint.dat");

上述语句可以把仿真到当前步的结果存储下来，然后利用loadBinaryBlock方法下面的语句进行重启动：

1
loadBinaryBlock(lattice, "checkpoint.dat");

使用此过程，一次只能写入一个数据块，即一次只能写入一个lattice。如果模拟的状态由各种块表示（例如，在使用几个块格的多组分流体中），则每个块都需要保存在单独的文件中。目前无法保存结构数据（Dynamic和数据处理器）。在加载模拟状态时，必须在加载数据之前重新创建上下文，例如边界条件。所以，这个*.dat文件其实就是一个仿真的初始值，只不过这个初始值来自于另一个仿真。

6 输入输出到文件

6.1 输出到终端和文件

在Palabos中，永远不要使用C++中常用的输入输出流对象cout，cerr或clog，或者更糟糕的是，使用printf将任何内容打印到终端。原因是这些方法都不是并行化的，如果有多个内核分别并行运行，会造成同时打印很多消息，造成死机。

这些符号由相应的可并行化版本pcout，pcerr和pclog替换。它们具有与传统语言相同的语法和相同的行为，但另外，它们在并行程序中提供了预期的行为：如果从一百个内核上并行化的程序打印消息，则消息只打印一次，而不是一百次。

可以使用这些并行化版本的输入输出流来打印字符串，数字，或者是提取自lattice的标量场，张量场，比如下面的语句：

1
2
3
4
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;

写入到文件：

1
2
plb_ofstream ofile("energy.dat");
ofile << setprecision(10) << *computeEnergy(lattice) << endl;

注意到，这里用的是plb_ofstream而不是标准的ofstream。同样的道理。

dat文件可以在MATLAB、Octave中使用。

文档中给了一个用Octave来处理的例子：

1
2
3
4
5
6
7
% Analysis of the data in Octave
% Assign manually the dimensions:
nx = 100;
ny = 100;
load energy.dat;
energy = reshape(energy, nx, ny);
imagesc(energy);

6.2 从文件中读取大数据集

通过类型为plb_ifstream的对象来读取来自于文件的数值，例如：

1
2
3
4
5
6
7
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;
}

注意，在使用MultiScalarField2D energy(nx,ny)的时候，要注意参数nx,ny和dat文件的维数相符。

但是，只有将数据流式传输到Palabos对象时，此方法才有效。换句话说，如果是你自己定义的C++对象，则不能使用plb_ifstream来操作。因为它会导致并行程序中出现意外行为。

解决这个问题有两个方法：

  1. 创建XML格式的参数文件，然后使用Palabos的自动XML解析工具对其进行解析；
  2. 仍然使用plb_ifstream读入文件，但是需要手动将数据传播到所有处理器，以便获得与数据并行兼容的工作并行程序Palabos的编程概念

6.2.1 采用XML格式的参数文件进行输入

文档中是推荐这个方法的，因为XML格式的参数文件结构良好，易读好懂。

Palabos提供了一种从XML文件中读取结构化输入数据的方法。此功能适用于串行和并行程序。然而，数据存储在普通的非并行数据结构中，这些数据结构在并行程序的所有线程上都是重复的。因此，XML文件应仅用于小型数据集，例如设置模拟所需的参数。，大数据集采用plb_ifstream读入。

Palabos能读入的XML文件有三个限制：

  1. 无法解析文档类型定义（DTD）或可扩展样式表语言（XSL）。
  2. 无法识别属性（但即使存在属性，仍然可以正确解析标记）。例如，在Palabos中无法访问以下XML选项卡中的属性id的值：。
  3. 给定的标签名称只能使用一次。如果多次使用，则只能在Palabos中访问最后一次出现。

这三个限制对不懂HTML语法的人来说，比如我，是没什么概念的。

但是可以依葫芦画瓢：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28



<Geometry>
    <inputFile> /home/user/data/inputFile.dat inputFile>
    <size>
    
        <nx> 10 nx>
        <ny> 20 ny>
        <nz> 30 nz>
    size>
    <viscosity>
    
        0.023
    viscosity>
    <umax>
    
        0.01
    umax>
Geometry>

<Inlet>
    
    <kind> Dirichlet kind>
    
    <values> 0 0.1 0.2 0.3 0.4 0.4 0.3 0.2 0.1 0 values>
Inlet>

在Palabos里面这样进行解析：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
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。

6.2.2 采用plb_ifstream

还是采用之前提到的energy.dat文件：

1
2
3
4
5
6
7
8
9
10
11
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;
}

注意到语句：global::mpi().bCast(&nx, 1);和global::mpi().bCast(&ny, 1);，这两个相当于是对全局进行了广播。这些个参数都被告知到了每个处理核中。这种处理输入值的方法在代码示例/ codesByTopic / io / manualUserInput.cpp中说明。

注意：

  • Palabos中没有用于交互式输入的对象xcin ; 用户输入必须通过输入文件或直接从命令行接收。

7 吐个槽，做个计划

暂时写到这里。大部分是翻译官方文档，没太多自己的东西。不过尽管如此，user guide搞了这么久还没搞定，需要加快速度。

搞完user guide之后，后面还要写一写边界条件的相关的类。这部分比较蛋疼，但是不着急。

数据处理器这部分也云里雾里，后续再看看。

再然后就是要看看Guo外力的动态类、以及Grid refinement。

实施IB-LBM在二维圆柱绕流上，这部分工作需要尽快完成。

时间比较紧，加油。

任务 所需时间
User Guide 2days(优先1)
进一步研究边界条件 2days(暂缓)
数据处理器 1days(优先1)
Grid refinement案例解读 1day(暂缓)
案例1 1day (优先2)
案例2 1day (优先2)
Guo类解读 1day(优先3)
圆柱绕流（IB-LBM） 5days(优先4)

严格执行，切勿纠缠！