Palabos User Guide中文解读 | 第十一章 | 输入/输出

作者的话：本人在学习palabos时，发现国内中文资料甚少，恰好网上可以直接搜到palabos user guide这种英文资料，加之时间充裕，便打算开始翻译，翻了一节后发现这可能算侵权，就比较伤脑筋，突然想到自己写中文解读即可，便有了下面的博客。

Palabos User Guide
Release 1.0 r1
Copyright © 2019 University of Geneva
Jul 05, 2019

Chapter Eleven

输入/输出

11.1 输出流动：编写至终端且写入文件

不要在palabos里用cout, ceer, clog, 或更糟糕，使用printf来输出信息到终端。它们都有对应的替代代码pcout，pceer，pclog，相同的语法和作用。更重要的是，在并行情况下也有理想的结果。

你可以使用它们来输出字符串、数字或者从lattice，标量场，张量场里抽取数据输出：

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替代：

plb_ofstream ofile(“energy.dat”);
ofile << setprecision(10) << *computeEnergy(lattice) << endl;

在交互式数据分析中（与matlab或Octave），向palabos里导入数据也很有用。

% Analysis of the data in Octave
% Assign manually the dimensions:
nx = 100;
ny = 100;
load energy.dat;
energy = reshape(energy, nx, ny);
imagesc(energy);

11.2 输入流动：从文件中读取大型数据集

输入文件和上节中谈到的输出文件差不多，通过plb_ifstrams类型的对象。
上节中提到的energy.dat可以通过下面的代码转化为矩阵的形式：

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格式文件，在11.4节 从XML文件读取用户输入 中可以找到解释。这样做不仅易于阅读，自动转换和纠错。

第二种方法，就是你不想使用XML格式文件时，你仍然可以用plb_ifstream，但你需要手动把数据添加到所有处理器中，这样确保并行运行，也符合palabos的数据并行程序概念。
我们假设有一个energy.dat文件，里面是nx*ny的矩阵，下面我们要导入它到palabos里，并自动生成尺寸正确的标量场：

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可以看到这种方法输入数据的代码。虽然可并行运算的程序也可单线程运行，你最好都这么做，确保两种方式运行不需要改程序。

在palabos里，从没有什么pcin可以用来输入，输入只有两种方式，要么文件要么代码。

11.3 访问代码行参数

和一般的C++程序一样，你可以通过main函数的argc和argv参数来访问代码行参数。

int main(int argc, char* argv[])

并行和单核运行效果一样。palabos里，这两个全局功能让你在程序任何位置都可以访问代码行参数。所以推荐你使用全局对象而不是main参数，它会自动转换类型和处理错误。如下面代码所示：

int main(int argc, char* argv[]) {
2 // Never forget to initialize Palabos, in every program.
3 plbInit(&argc, &argv);
45 pcout << "The number of input arguments is: " 6 << global::argc() << endl;
78 double double_argument;
9 plint plint_argument;
10 string string_argument;
11
12 try {
13 global::argv(1).read(double_argument);
14 global::argv(2).read(plint_argument);
15 global::argv(3).read(string_argument);
16 }
17 // React to errors, either if there are less than 3 arguments,
18 // or if one of the arguments fails to convert to the expected
19 // data type.
20 catch(PlbIOException& exception) {
21 // Print the corresponding error message.
22 pcout << exception.what() << endl;
23 // Terminate the program.
24 exit(1);
25 }
26
27 // … Execute the main program …
28 }

在这里，有个C++的异常机制被用于处理用户的错误输入。假若你不想处理这些错误，你可以用read函数的noThrow版本。

bool ok = global::argv(1).readNoThrow(double_argument);
if (!ok) {
pcout << “Error while reading argument 1.” << endl;
}

例子：examples/showCases/rayleighTaylor2D.cpp和rayleighTaylor3D.cpp。

11.4 从XML文件读取用户输入

并行和单核运行情况下，palabos都可以读取XML文件的输入数据。数据会在每个线路都复制一份。

由此XML文件一般都是小的数据集，如初始化状态的参数设置。

大的数据集，如lattice，应该直接被读入palabos的数据结构中，从而数据被分配到并行机制中，如11.2节 输入流动：从文件中读取大型数据集 所示。

开源库TinyXML可以帮助你处理语法结构，这个库无需你安装，palabos软件包自带。
虽然palabos可以读取任何正常生成的XML文件，但也有3个限制：
1）DTDs或XSL不可行。

2） Attributes are not recognized (but a tag is still properly parsed, even if attributes are present). For example, thebvalue of the attribute id in the following XML tab cannot be accessed in Palabos: < someTag id=“5”>.

3）一个标签名只能使用一次，多次使用以最后一次的数据为准。

下面是一个palabos可用的典型的输入文件：
XML file myInput.xml:

1<?xml version="1.0" ?>
2 < !-- Configuration parameters for my simulation --> 3
4 < Geometry>
5 < inputFile> /home/user/data/inputFile.dat < /inputFile>
6 < size>
7 < !-- Number of lattice nodes in each space direction. -->
8 < nx> 10 < /nx> 9 < ny> 20 < /ny>
10 < nz> 30 < /nz>
11 < /size>
12 < viscosity>
13 < !-- Viscosity in lattice units. -->
14 0.023
15 < /viscosity>
16 < umax>
17 < !-- Maximum velocity in lattice units. -->
18 0.01
19 < /umax>
20 < /Geometry>
21
22 < Inlet>
23 < !-- Use a velocity Dirichlet boundary condition. -->
24 < kind> Dirichlet < /kind>
25 < !-- Velocity profile values on the inlet, in
26 dimensionless units. -->
27 < values> 0 0.1 0.2 0.3 0.4 0.4 0.3 0.2 0.1 0 < /values>
28 < /Inlet>

下面是对应的palabos代码：

1 try {
2 // Open the XMl file.
3 XMLreader xmlFile(“myInput.xml”);
4
5 // It’s good policy to flush the content of the XML
6 // file to the terminal, so that in future, when
7 // you check the program’s output, you remember
8 // which parameters where used.
9 pcout << “CONFIGURATION” << endl;
10 pcout << “=============” << endl << endl;
11 xmlFile.print(0);
12
13 string inputFile;
14 xmlFile[“Geometry”][“inputFile”].read(inputFile);
15
16 plint nx, ny, nz;
17 xmlFile[“Geometry”][“size”][“nx”].read(nx);
18 xmlFile[“Geometry”][“size”][“ny”].read(ny);
19 xmlFile[“Geometry”][“size”][“nz”].read(nz);
20
21 T viscosity, uMax;
22 xmlFile[“Geometry”][“viscosity”].read(viscosity);
23 xmlFile[“Geometry”][“umax”].read(uMax);
24
25 string inletKind;
26 xmlFile[“Inlet”][“kind”].read(inletKind);
27 vector inletValues;
28 xmlFile[“Inlet”][“values”].read(inletValues);
29 }
30 catch (PlbIOException& exception) {
31 pcout << exception.what() << endl;
32 return -1;
33 }

从一个标签名中你可以读取整个数据行列，就像上例最后一个标签名所示的一样，会被录入向量 inletValues。

11.5 生成2D和3D图像

你完全可以生成整个流域，或者小区域，或者标量场的一个切片的图像。默认情况下palabos输入PPM格式的图像，一种无需其他外部图像库帮助即可生成的ASCII格式。如果你安装了ImageMagick包（代码行的工具covert就在这个包里），就可以输出gif格式。

第一步，你需要创建一个图像输出的对象：

ImageWriter imageWriter(“leeloo”);

leeloo代表着标量变量转为RGB色彩模式的colormap。除了leeloo还有earth，water，air，fire。

下一步你可以看到各种选项：

// 在[minValue, maxValue]范围内输出PPM格式的图像.
imageWriter.writePpm(scalarField, “myImage”, minValue, maxValue);

// 在自动生成的范围内输出PPM格式的图像.
imageWriter.writeScaledPpm(scalarField, “myImage”);

// 在[minValue, maxValue]范围内输出GIF格式的图像.
imageWriter.writeGif(scalarField, “myImage”, minValue, maxValue);

// 在自动生成的尺寸范围内输出PPM格式的图像
imageWriter.writeScaledGif(scalarField, “myImage”);

// 在自动生成尺寸的范围内输出GIF格式的图像，并重新调整尺寸至x*y。（纵横比不变）
imageWriter.writeScaledGif(scalarField, “myImage”, sizeX, sizeY);

这种方法也可以生成3D模拟的截图。
即切片一下3D模拟的数据，即选择一个仅有一格子高的维度。

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

11.6 VTK输出的后处理

任何关于block-lattice，标量场，张量场的数据都可以被VTK格式导出。此时你需要了解paraview。

尽管后处理不需要跑程序时那么要求精确度，但使用paraview的确是非常精确。若是你想要节省一些内存，这时可以把双精度转化为单精度，正如下面这个3D例子所示：

// 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< T> 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文件时，你想添加多少block都可以。多少个输出参数都是自定义的，上面的例子里有3个张量场，1个标量场。

2D数据是用VtkImageOutput2D来输出的，例子的路径：examples/showCases/poiseuille/。

paraview处理2D数据的分析不是最好的选择，你可以通过11.1节 _输出流动：编写至终端且写入文件 _ 将数据写入文本文件，再用Octave或matlab（Python、R、SciLab/ OO-Spreadsheet等一切你喜欢的）来实现可视化。

11.7 检查点：保存和载入模拟情况

尽管你可以使用对象VtkImageOutputXD来写出VTK文件，但你无法把VTK文件的数据反向导入palabos里。

如果你希望保存运行状态，以便于未来恢复，可以使用函数saveBinaryBlock，可在例子examples/codesByTopic/io/checkPointing.cpp中看到展示。

saveBinaryBlock(lattice, “checkpoint.dat”);

随后可以使用loadBinaryBlock函数载入:

loadBinaryBlock(lattice, “checkpoint.dat”);

这种情况下，你只需要一个block写入一次数据。如果模拟是多个不同block一起运行的（比如一些block-lattice一起运行的多组流），每个block都需要分别保存一下。

目前不可以保存结构性数据（dynamics对象和数据处理器）。当载入模拟状态时，一些内容，比如边界条件，需要在导入数据之前重新生成。