Palabos用户手册翻译及学习（三）设置中的问题与定义边界条件

设置中的问题与定义边界条件

Palabos用户文档 的第八章和第九章

（一）Setting Up A Problem

原始文档

Attributing dynamics objects

When constructing a new block-lattice, you must decide what type of collision is going to be executed on each of its cells, by assigning them a dynamics object. To avoid bugs related to cells without dynamics objects, the constructor of a block-lattice assigns a default-alue for the dynamics to all cells, the so-called background dynamics:

// Construct a new block-lattice with a background-dynamics of type BGK.
MultiBlockLattice3D<T,DESCRIPTOR> lattice( nx, ny, nz, new BGKdynamics<T,DESCRIPTOR>(omega) );

After this, the dynamics of each cell can be redefined in order to adjust the behavior of the cells locally. The background dynamics differs from usual dynamics objects in an important way. To obtain memory savings, the constructor of the block-lattice creates only one instance of the dynamics object, and all cells refer to the same instance. If you modify a dynamics object on one cell, it changes everywhere. As an example, consider the relaxation parameter omega, which is stored in the dynamics, not in the cell. A function call like

lattice.get(0,0,0).getDynamics().setOmega(newOmega); 

would affect the value of the relaxation parameter on each cell of the lattice. Note that this particular behavior of the background dynamics, having several cells referring to the same object, cannot be reproduced by other means than constructing a new block-lattice. All other functions which override the background-dynamics on given cells with a new dynamics object are defined in such a way as to create an independent copy of the dynamics object for each cell.
If the reference semantics of the background dynamics is contrary to your needs, you can re-define all dynamics objects right after constructing a block-lattice, to guarantee that all cells have an independent copy:

// Override background-dynamics to guarantee and independent per-cell
// copy of the dynamics object.
defineDynamics( lattice, lattice.getBoundingBox(), new BGKdynamics<T,DESCRIPTOR> );

There exist different variants of the function defineDynamics which you can use to adjust the dynamics of a group of cells (their syntax can be found in the Appendix Operations on the block-lattice):

• One-cell version: Assign a new dynamics to just one cell.
  Example: tutorial/tutorial2/tutorial2_3.cpp.
• BoxXD version: Assign a new dynamics to all cells within a rectangular domain.
  Example: showCases/multiComponent2d/rayleighTaylor2D.cpp, or the 3d example.
• DotListXD version: Assign a new dynamics to several cells, listed individually in a dot-list structure.
  Example: codesByTopic/dotList/cylinder2d.cpp
• Domain-functional version: Provide an analytical function which indicates the coordinates of cells which get a new dynamics object.
  Example: showCases/cylinder2d/cylinder2d.cpp
• Bool-mask version: Specify the location of cells which get a new dynamics object through a Boolean mask, represented through a scalar-field.
  Example: codesByTopic/io/loadGeometry.cpp
  
  

Initial values of density and velocity

It is quite common to assign an initial value of velocity and density to a lattice by initializing all cells to an equilibrium distribution with the chosen density and velocity value. While this approach is not sufficient for time-dependent benchmark problems which depend critically on the initial state, it is most often sufficient to get a simulation started with reasonable initial values. The function initializeAtEquilibrium (see Appendix Operations on the blocklattice) is provided to initialize the cells within a rectangular-shaped domain at an equilibrium distribution. It comes in two flavors: one with a constant density and velocity within the domain (see the example examples/showCases/cavity2d or 3d), and one with a space-dependent value of these macroscopic values, specified through a userdefined function (see the example examples/showCases/poiseuille).
To initialize cells in a different way, it is simplest to apply a custom operator to the cells with the help of one-cell functionals and indexed one-cell functionals introduced in Section Convenience wrappers for local operations.

文档翻译

设置动态对象

当构造一个新的块格时，你必须通过分配一个动态对象来决定在每个单元格上执行什么类型的碰撞。为了避免与没有动态对象的单元格相关的bug，块格的构造器将动态的默认值赋给所有单元格，即所谓的后台动态:

// 以BGK动力学类作为背景动力学构造一个新的块格。
MultiBlockLattice3D<T,DESCRIPTOR> lattice( nx, ny, nz, new BGKdynamics<T,DESCRIPTOR>(omega) );

在这之后，每个单元格的动态可以重新定义，以调整单元格的局部行为。背景动力学与一般的动力学对象有重要的区别。为了节省内存，块格的构造函数只创建一个动态类对象的实例，所有单元格都引用同一个实例。如果您在任一个单元格上修改动态类对象，它将到处更改。例如，考虑松弛参数，它存储在动态类对象中，而不是在单元格中。如以下函数调用

lattice.get(0,0,0).getDynamics().setOmega(newOmega); 

会影响晶格中每个单元的弛豫参数的值。请注意，这是背景动力学的特定行为，几个单元格指向同一个对象，不能通过其他方法重制，只能通过构造一个新的块晶格来实现。使用新的动态类对象覆盖给定单元上的背景动态的所有其他函数的定义方式是，为每个单元创建动态类对象的独立副本。
如果背景动态的引用语义与你的需要相反，你可以在构建一个块格之后重新定义所有的动态对象，以保证所有的单元都有一个独立的副本:

// 覆盖背景动态以保证每个单元都有动态对象的独立副本。
defineDynamics( lattice, lattice.getBoundingBox(), new BGKdynamics<T,DESCRIPTOR> );

你可以使用defineDynamics函数的不同变体来调整一组细胞的动态(它们的语法可以在block-lattice的附录操作中找到):

• One-cell版本：将一个新的动态分配给一个细胞。
  例子：tutorial/tutorial2/tutorial2_3.cpp
• BoxXD版本：将一个新的动态分配给矩形域中的所有单元格。
  例子：showCases/multiComponent2d/rayleighTaylor2D.cpp或 3d 案例
• 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)。

解释说明

在构建块格时生成的动态类对象只有一个，在堆区分配对象空间，然后给块格中的每一个单元传递一个动态类对象的指针，所以只有一个对象，所有单元格共享。为特定单元格单独设置动态类，是通过在兑取新建一个相应的动态类对象并将其指针传递给对应单元格来覆盖只想背景动态类对象的指针，从而达到重新设置类对象的效果。One-cell、BoxXD、DotListXD、Domain-functional、Bool-mask都是指定更改所限定作用范围的参数，你可以把它们都想象成stl标准模板库中的iterator迭代器，不仅在动态类的设定中有使用，基本上所有其他的设置，都有他们的身影。
密度和速度设置，整个块格的密度和速度的初始化，方腔流是将整个流区初始化为相同的密度和速度，也就是稳态模拟最有可能达到的平衡态附近，同理，泊肃叶流是将整个流区初始化为泊肃叶分布形式的密度和速度，也是对应的稳态模拟最可能达到的平衡态的附近，此种设置，都是为了能让模拟更快达到稳态结果，如果不这样设置，可能会导致模拟失败，非稳态的计算，我认为这个初始值并不是特别重要，主要的是在边界以及特定流动区域不要出现计算的死点，才是初始值和对应模拟参数需要关注的重点。

（二）Defining Boundary Condition

原始文档

Overview

The library Palabos currently implements a bunch of different boundary conditions for straight, grid-aligned boundaries. For general boundaries, the available options include stair-cased walls through bounce-back nodes, and higherorder accurate curved boundaries.
IMPORTANT: Curved boundaries are available as of version 1.0, and are part of a fully-featured parallel stack, including parallel voxelization and I/O extensions. Curved boundaries are not yet documented. You can however find a full-featured example in showCases/aneurysm.

Grid-aligned boundaries

The class OnLatticeBoundaryConditionXD
Some of the implemented boundary conditions are local (and are therefore implemented as dynamics classes), some are non-local (and are therefore implemented as data processors), and some have both local and nonlocal components. To offer a uniform interface to the various possibilities, Palabos offers the interface OnLatticeBoundaryConditionXD which is responsible for instantiating dynamics objects, adding data processors, or both, depending on the chosen boundary condition. The following line for example is used to chose a Zou/He boundary condition:

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createZouHeBoundaryCondition3D<T,DESCRIPTOR>();

The four possibilities currently offered are (see palabos.org-bc for details):

model	class
Regularized BC	createLocalBoundaryConditionXD
Skordos BC	createInterpBoundaryConditionXD
Zou/He BC	createZouHeBoundaryConditionXD
Inamuro BC	createInamuroBoundaryConditionXD

With each of these boundary conditions, it is possible to implement velocity Dirichlet boundaries (all velocity components have an imposed value) and pressure boundaries (the pressure is imposed, and the tangential velocity components are zero). Furthermore, various types of Neumann boundary conditions are proposed. In these case, a zero-gradient for a given variable is imposed with first-order accuracy by copying the value of this variable from a neighboring cell.

Boundaries on a regular block domain
The usage of these boundary conditions is somewhat tricky, because Palabos needs to know if a given wall is a straight wall or a corner (2D), or a flat wall, an edge or a corner (3D), and it needs to know the wall orientation. To simplify this, all OnLatticeBoundaryConditionXD objects have a method setVelocityConditionOnBlockBoundaries and a method setPressureConditionOnBlockBoundaries, to set up boundaries which are located on a rectangular domain. If all nodes on the boundary of a given 2D box boundaryBox implement a Dirichlet condition, use the following command:

Box2D boundaryBox(x0,x1, y0,y1);
boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, boundaryBox, locationOfBoundaryNodes );

The first argument indicates the boxed shape on which the boundary is located, and the second argument specifies on which node a velocity condition is to be instantiated. If all boundary nodes are located on the exterior bounding box of the lattice, the above command simply becomes:

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);

Finally, if simply all nodes of the exterior bounding box shall implement a velocity condition, this simplifies to:

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);

Now, suppose that the upper and lower walls, including corners, in a nx -by- ny lattice implement a free-slip condition (zero-gradient for tangential velocity components), the left wall a Dirichlet condition, and the right wall and outflow condition (zero-gradient for all velocity components). This is achieved by adding one more parameter to the function call:

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 );

Remember that if the boundary is not located on the outer bounds of the lattice, you must use an additional Box3D argument to say where the boundaries are located, additionally to saying which boundary is going to be added. This is necessary, because Palabos needs to know the orientation and the nature (flat wall, corner, etc.) of the boundary nodes.

Zero-gradient conditions
The optional arguments like boundary::outflow are of type boundary::BcType, and can have the following values for velocity boundaries:

object values	explain
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

For pressure boundaries, you have the following choice:

object values	explain
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

Boundary conditions cannot be overriden
It is not possible to override the type of a boundary. Once a boundary node is a Dirichlet velocity node, it stays a Dirichlet velocity node. The result of re-defining it as, say, a pressure boundary node, is undefined. Therefore, boundary conditions must be defined carefully, and piece-wise, as shown in the example above. On the other hand, the value imposed on the boundary (i.e. the velocity value on a Dirichlet velocity boundary) can be changed as often as needed.

Setting the velocity or pressure value on Dirichlet boundaries
A constant velocity value is imposed with the function setBoundaryVelocity. The following line sets the velocity values to zero on all nodes of a 2D domain which have previously been defined to be Dirichlet velocity nodes:

setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );

On all other nodes, this command has no effect. To set a constant pressure value (or equivalently, density value), use the command setBoundaryDensity:

setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. );

A non-constant, space dependent velocity resp. pressure profile can be defined by replacing the velocity resp. pressure argument by either a function or by a function object (an instance of a class which overrides the function call operator) which yields a value of the velocity resp. density as a function of space position. An example is provided in the file examples/showCases/poiseuille/poiseuille.cpp.

Interior and exterior boundaries
Sometimes, the boundary is not located on the surface of a regular, block-shaped domain. In this case, the functions like setVelocityConditionOnBlockBoundaries are useless, and the boundary must be manually constructed. As an example, here’s how to construct manually a free-slip condition along the top-wall, including the corner nodes:

// 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 );

A distinction is made between external and internal corners, depending on whether the corner is convex or concave. On the following example geometry, you’ll find five external and one internal corner:

The extensions like 1P, NP, and PP at the end of the methods of the boundary-condition object are used to indicate the orientation of the wall normal, pointing outside the fluid domain, as shown on the figure above. On a straight wall, the code 1P means: “the wall normal points into positive y-direction”. Likewise, the inlet would be labeled with the code 0N as in “negative x-direction”. On a corner, the code NP means “negative x-direction and positive y-direction”. It is important to mention that this wall normal is purely geometrical and does not depend on whether a given wall has the function of an inlet or an outlet. In both cases, the wall normal points away from the fluid, into the non-fluid area.
In 3D, you’ll use the following function for plane walls:

addVelocityBoundaryDO

where the direction D can be 0 for x, 1 for y, or 2 for z, and the orientation O has the value P for positive and N for negative. Edges can be internal or external. For example:

addInternalVelocityEdge0NP
addExternalVelocityEdge1PN

where the first digit of the code indicates the axis to which the edge is parallel, and the two subsequent digits indicate the orientation of the edge inside the plane normal to the edge, in the same way as the corner nodes in 2D. The axes are counted periodically: 0NP means x-plane, negative y-value, and positive z-value, whereas 1PN means y-plane, positive z-value, and negative x-value. Finally, 3D corners are constructed in the same way as in 2D, through a function call like the following:

addInternalVelocityCornerPNP
addExternalVelocityCornerNNN

Periodic boundaries

The concept of periodicity in Palabos is different from the approach chosen for the other types of boundary conditions. Periodicity can only be implemented on opposite, outer boundaries of an atomic-block or a multi-block. This behavior works also if the multi-block has a sparse memory implementation. In this case, all fluid nodes which are in contact with an outer boundary of the multi-block can be periodic, if the corresponding node on the opposite wall is a fluid node.
In the case of an atomic-block (remember that this case is uninteresting, because you should always work with multiblocks anyway), periodicity is simply a property of the streaming operator. It has no effect for scalar-fields and tensor-fields. In the case of a multi-block however, periodicity can also affect scalar-fields and tensor-fields. Being periodic in this case means that if you define a non-local data processor which accesses nearest neighbor nodes, the value of these neighbor nodes is determined by periodicity along an outer boundary of the multi-block.
By default, all blocks are non-periodic: as mentioned previously, the behavior of outer boundary nodes defaults to one of the versions of bounce-back algorithm. To turn on periodicity along, say, the x-direction, write a command like

lattice.periodicity().toggle(0, true); // Periodic block-lattice.
scalarField.periodicity().toggle(0, true); // Periodic scalar-field.

You can also define all boundaries to be periodic through a single command:

lattice.periodicity().toggleAll(true);

Please note that the periodicity or non-periodicity of a block should be defined as soon as possible after the creation of the block, because Palabos needs to know if boundaries are periodic or not when adding or executing data processors.
As a last remark, you should be aware that boundaries of Neumann type (outflow, free-slip, adiabatic thermal wall, etc.) should never be defined to be periodic. This wouldn’t make sense anyway, and it leads to an undefined behavior (aka program crash).

Bounce-back

In Palabos, all outer boundaries of a lattice which are not periodic automatically implement a version of the bounceback boundary algorithm, according to the following algorithm. In pre-collision state, all unknown populations are assigned the post-collision value on the same node at the end of the previous iteration, from the opposite location. We will call this algorithm Bounce-Back 1.
A bounce-back condition can also be used elsewhere in the domain, by explicitly assigning a dynamics object of type BounceBack<T,Descriptor> to chosen cells. On these newly assigned nodes, the collision step is replaced by a bounce-back procedure, during which each population is replaced by the population corresponding to the opposite direction. Such a node can not be considered any more as a fluid node. Computing velocity-moments of the populations on a bounce-back node yields for example arbitrary results, because some of the populations (those which are not incoming from the fluid) have arbitrary values. Consequently, you cannot compute macroscopic quantities like density and velocity on these nodes. Instead, these nodes should be considered as pure algorithmic entities which have the purpose to re-inject into the fluid all populations that leave the domain. Consider fluid nodes adjacent to this type of bounce-back cells. If you think about it, you’ll realize that these nodes behave exactly as if they were cells of type Bounce-Back 1, with a small difference: unknown populations at time t get the post-collision value from opposite populations from time t-2 (and not time t-1 as in Bounce-Back 1. We will refer to this version of bounce-back as Bounce-Back 2.
The effect of both versions of bounce-back is to produce a no-slip wall located half-way between nodes. In the case of Bounce-Back 1, this is half a cell-width beyond the bounce-back fluid cell, whereas in the case of Bounce-Back 2 it is half-way between the last fluid cell and the non-fluid bounce-back node. Obviously, Bounce-Back 2 wastes numerical precision over Bounce-Back 1 by leaping over a time step, and therefore leads to a lower-order accuracy in time of the numerical method. Furthermore, both versions of bounce-back represent the wall location through a staircase approximation, which is low-order accurate (“Order-h”) in space. On the other hand, Bounce-Back 2 offers a huge advantage: its implementation is independent of the wall orientation. If you decide that a node is a bounce-back node, you simply say so; no need to know if it is a corner, a left wall, or a right wall. Constructing a geometry is then as easy as filling areas with bounce-back nodes. In many cases, this ease of use makes up, to a large extent, for the low-order accuracy of the wall representation. When setting up a new simulation, it is always good to try Bounce-Back 2 as a first attempt, as the quality of the obtained results is often surprising. In highly complex geometry such as porous media, bounce-back is often even considered superior to other approaches.

Bounce-back domain from analytical description
The Rule 0 in Palabos can be stated often enough: don’t write loops over space in end-user code. This also true, of course, when assigning bounce-back dynamics to lattice cells. To guarantee reasonable performance, this is task is performed inside a data processor, or even better, by using the wrapper function defineDynamics. This process is illustrated in the example program examples/codesByTopic/bounceBack/instantiateCylinder.cpp. First, a class is written which defines the location of the bounce-back nodes as a function of the cell coordinates iX and iY. In the present case, the nodes are located inside a 2D circular domain:

template<typename T>
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<T>* clone() const {
		return new CylinderShapeDomain2D<T>(*this);
	}
private:
	plint cx, cy;
	plint radiusSqr;
};

An instance of this class is then provided as an argument to defineDynamics:

defineDynamics( lattice, lattice.getBoundingBox(),
				new CylinderShapeDomain2D<T>(cx, cy, radius),
				new BounceBack<T,DESCRIPTOR> );

The second argument is of type Box2D, and it can be used to improve the efficiency of this function call: the bounce-back nodes are attributed on cells on the intersection between this box and the domain specified by the domain-functional.
If your domain is specified as the union of a collection of simpler domains, you can use defineDynamics iteratively on each of the simpler domains. If the domain is the intersection or any other geometric operation of simpler domains, you’ll need to play with boolean operators inside the definition of your domain-functional.

Bounce-back domain from boolean mask
In this section, it is assumed that the geometry of the domain is prescribed from an external source. This is for example the case when you simulate a porous media and possess data from a scan of the porous material in question. This easiest approach consists in storing this data as an ASCII file which distinguishes fluid nodes from solid nodes by zeros and ones, separated with spaces, tabs, or newline characters. The file represents the content of a regular array and stores only raw data, no information on the matrix dimensions. The data layout must be the same as in Palabos: an increment of the z-coordinate represents a continuous progress in the file (in 3D), followed by the y-coordinate, and then the x-coordinate. This data is simply flushed from the file into a scalar-field. The following code is taken from the example examples/codesByTopic/io/loadGeometry.cpp:

MultiScalarField2D<T> boolMask(nx, ny);
plb_ifstream ifile("geometry.dat");
ifile >> boolMask;

Note that currently, no error-checking is implemented for such I/O operations. It is your responsibility to ensure that the dimensions of the scalar-field correspond to the size of the data in the geometry file. As a next step, the bounce-back nodes are instantiated using one of the versions of the function defineDynamics:

defineDynamics(lattice, boolMask, new BounceBack<T,DESCRIPTOR>, true);

This function is currently only defined for multi-blocks using the same data type (T in this case), which explains why the scalar-field boolMask is of type T and not bool, as one could have expected. The scalar-field is not bool, as one could have expected, but the real type T.

文档翻译

综述

Palabos库目前为网格对齐的直线边界实现了许多不同的边界条件。对于一般的边界，可用的选项包括通过反弹节点的阶梯式墙壁，以及更高阶的精确弯曲边界。
重要提示:从1.0版开始就可以使用曲线边界，它是一个功能齐全的并行堆栈的一部分，包括并行体素化和I/O扩展。弯曲的边界还没有记录。然而，你可以在showCases/aneurysm中找到一个全功能的例子。

网格对齐边界

OnLatticeBoundaryConditionXD类
一些实现的边界条件是本地的(因此被实现为动态类)，一些是非本地的(因此被实现为数据处理器)，一些同时具有本地和非本地组件。为了提供各种可能性的统一接口，Palabos提供OnLatticeBoundaryConditionXD接口，该接口负责实例化动态对象、添加数据处理器，或两者兼得，这取决于所选择的边界条件。以下面这行为例，选择一个邹/河边界条件:

OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition = createZouHeBoundaryCondition3D<T,DESCRIPTOR>();

目前提供的四种可能性是:

model	class
Regularized BC	createLocalBoundaryConditionXD
Skordos BC	createInterpBoundaryConditionXD
Zou/He BC	createZouHeBoundaryConditionXD
Inamuro BC	createInamuroBoundaryConditionXD

对于每一个边界条件，都可以实现速度Dirichlet边界(所有速度分量都有一个施加的值)和压力边界(施加压力，切向速度分量为零)。此外，还提出了各种类型的Neumann边界条件。在这种情况下，通过从邻近单元格复制该变量的值，以一阶精度对给定变量施加零梯度。

常规块域上的边界
这些边界条件的使用有些棘手，因为Palabos需要知道给定的墙是直墙还是角(2D)，还是平墙、边或角(3D)，并且需要知道墙的方向。为了简化这一点，所有OnLatticeBoundaryConditionXD对象都有一个方法setVelocityConditionOnBlockBoundaries和一个方法setPressureConditionOnBlockBoundaries，用来设置位于矩形域中的边界。如果给定二维框边界上的所有节点都实现Dirichlet条件，则使用以下命令:

Box2D boundaryBox(x0,x1, y0,y1);
boundaryCondition->setVelocityConditionOnBlockBoundaries ( lattice, boundaryBox, locationOfBoundaryNodes );

第一个参数表示边界所在的方箱形状，第二个参数指定要在哪个节点上实例化velocity条件。如果所有边界节点都位于格外边界框内，则上述命令简单为:

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice, locationOfBoundaryNodes);

最后，如果外边界框的所有节点都要实现速度条件，则简化为:

boundaryCondition->setVelocityConditionOnBlockBoundaries(lattice);

现在，假设在一个二维晶格中，上下壁面(包括角)实现了一个自由滑动条件(切向速度分量为零梯度)，左壁面为一个Dirichlet条件，右壁面和出流条件(所有速度分量为零梯度)。这是通过在函数调用中增加一个参数来实现的:

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需要知道边界节点的方向和性质(平墙、角等)。

零梯度条件
The optional arguments like boundary::outflow are of type boundary::BcType, and can have the following values for velocity boundaries:
可选参数如boundary::outflow的类型为boundary::BcType，速度边界可以有以下值:

object values	explain
boundary::dirichlet (default value)	Dirichlet 条件： 实施速度值
boundary::outflow or boundary::neumann	所有速度分量零梯度
boundary::freeslip	切向速度分量为零梯度，其它分量为零
boundary::normalOutflow	法向速度分量为零梯度，其他分量为零

对于压力边界,你有以下选择:

object values	explain
boundary::dirichlet (default value)	Dirichlet 条件： 施加压力值，切向速度分量为零
boundary::neumann	压强为零梯度，切向速度分量为零

边界条件不能被覆盖
不可能重写边界的类型。一旦边界节点是狄利克雷速度节点，它仍然是狄利克雷速度节点。将其重新定义为，例如，一个压力边界节点的结果是未定义的。因此，边界条件必须谨慎地、分段地定义，如上例所示。另一方面，施加在边界上的值(即狄利克雷速度边界上的速度值)可以根据需要随时改变。

在狄利克雷边界上设置速度或压力值
由函数setBoundaryVelocity施加一个恒定的速度值。下面的线将二维域中所有节点上的速度值设为0，之前定义为狄利克雷速度节点:

setBoundaryVelocity(lattice, lattice.getBoundingBox(), Array<T,2>(0.,0.) );

在所有其他节点上，此命令无效。要设置一个恒定的压力值(或等效的密度值)，使用命令setBoundaryDensity:

setBoundaryDensity(lattice, lattice.getBoundingBox(), 1. );

A non-constant, space dependent velocity pressure profile can be defined by replacing the velocity resp. pressure argument by either a function or by a function object (an instance of a class which overrides the function call operator) which yields a value of the velocity resp. density as a function of space position. An example is provided in the file examples/showCases/poiseuille/poiseuille.cpp.
一个非恒定的，依赖于空间的速度分布和压力分布可以分别通过替换速度分布和压力分布来确定。一个函数或一个函数对象(一个覆盖函数调用操作符的类的实例)的压力参数，该参数分别产生速度和的值作为空间位置的函数。示例在examples/showCases/poiseuille/poiseuile .cpp文件中提供。

内外边界
有时，边界并不位于规则的块状区域的表面。在这种情况下，像setVelocityConditionOnBlockBoundaries这样的函数是没有用的，必须手工构造边界。作为一个例子，这里展示如何手动构建一个沿顶壁自由滑动的边界条件，包括角节点:

// 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 );

外部和内部的角是有区别的，这取决于角是凸的还是凹的。在下面的几何示例中，你会发现五个外部和一个内部角落:

如上图所示，边界条件对象方法末尾的1P、NP和PP等扩展用于指示壁面法线指向流体域外的方向。在一面直墙上，代码1P表示:“墙法线指向正的y方向”。同样地，入口将被标记代码0N为“负x方向”。在一个角落里，NP代表“负的x方向和正的y方向”。需要指出的是，这道墙法线完全是几何图形，并不取决于给定的墙是具有入口还是出口的功能。在这两种情况下，壁法线点远离流体，进入非流体区。
在3D中，你将使用以下功能来制作平面墙:

addVelocityBoundaryDO

方向D可以是x为0,y为1,z为2，方向O为P为正，N为负。边可以是内部的，也可以是外部的。例如:

addInternalVelocityEdge0NP
addExternalVelocityEdge1PN

其中，代码的第一个数字表示与边缘平行的轴，随后的两个数字表示与边缘垂直的平面内边缘的方向，与2D中的角节点相同。轴是周期性计数的:0NP表示x平面，负的y值，正的z值，而1PN表示y平面，正的z值，负的x值。最后，通过如下所示的函数调用，以与2D相同的方式构造3D角:

addInternalVelocityCornerPNP
addExternalVelocityCornerNNN

周期边界

Palabos中周期性的概念不同于其他类型边界条件所选择的方法。周期性只能在原子块或多块的相对、外部边界上实现。如果多块具有稀疏内存实现，则此行为也有效。在这种情况下，如果边界对面壁上的对应节点是流体节点，则与多块外边界接触的所有流体节点都是周期性的。
在原子块的情况下(请记住，这种情况是无趣的，因为无论如何都应该使用多块)，周期性只是流操作符的一个属性。它对标量场和张量场没有影响。然而，在多块的情况下，周期性也会影响标量域和张量域。在这种情况下，周期性意味着如果定义一个访问最近邻居节点的非本地数据处理器，这些邻居节点的值将由沿多块外部边界的周期性决定。
默认情况下，所有块都是非周期性的：如前所述，外部边界节点的行为默认为一个反弹算法的版本。要打开沿x方向的周期性，可以这样编写命令

lattice.periodicity().toggle(0, true); // Periodic block-lattice.
scalarField.periodicity().toggle(0, true); // Periodic scalar-field.

你也可以定义所有的边界是周期性的通过一个单一的命令:

lattice.periodicity().toggleAll(true);

请注意，块的周期性或非周期性应该在块创建后尽快定义，因为Palabos需要知道添加或执行数据处理器时边界是否是周期性的。
最后，你应该知道Neumann类型的边界(流出、自由滑动、绝热热壁等)不应该被定义为周期性的。无论如何，这都没有意义，并且会导致未定义的行为(即程序崩溃)。

反弹

在Palabos中，晶格中所有非周期的外边界都会根据下面的算法自动实现一个版本的反弹边界算法。在碰撞前状态下，所有未知种群在前一次迭代结束时，从相反的位置在同一节点上分配碰撞后值。我们称这个算法为Bounce-Back 1。
通过显式地将类型为BounceBack<T,Descriptor> 分配给所选的单元格，还可以在域中的其他地方使用反弹条件。在这些新分配的节点上，碰撞步骤将被一个反弹过程替换，在此过程中，每个种群将被对应于相反方向的种群替换。这样的节点不能再被认为是流动节点。例如，在回弹节点上计算种群的速度矩可以得到任意的结果，因为有些种群(那些不是来自流体的种群)具有任意的值。因此，你不能在这些节点上计算像密度和速度这样的宏观量。相反，这些节点应该被视为纯粹的算法实体，其目的是将所有离开域的种群重新注入流体。考虑与这种弹回单元相邻的流体节点。如果你仔细想想，你会发现这些节点的行为与Bounce-Back 1类型的单元完全一样，只是有一点不同：t时刻的未知种群从t-2时刻的相反种群获得碰撞后的值(而不是像弹回1那样从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中的规则0需要经常声明：不要在最终用户代码中编写空间上的循环。当然，当将反弹动力学分配给晶格单元时，也是如此。为了保证合理的性能，这是在数据处理器中执行的任务，或者更好，使用包装器函数defineDynamics。这个过程在示例程序examples/codesByTopic/bounceBack/instantiateCylinder.cpp中进行了说明。首先，编写一个类，该类将反弹节点的位置定义为单元坐标iX和iY的函数。在本例中，节点位于二维圆形域内:

template<typename T>
class CylinderShapeDomain2D : public DomainFunctional2D {
public:
	CylinderShapeDomain2D(plint cx_, plint cy_, plint radius) 
						: cx(cx_), cy(cy_), radiusSqr(util::sqr(radius)) { }
	// 函数调用操作重写被指定反弹节点的位置。
	virtual bool operator() (plint iX, plint iY) const {
		return util::sqr(iX-cx) + util::sqr(iY-cy) <= radiusSqr;
	}
	virtual CylinderShapeDomain2D<T>* clone() const {
		return new CylinderShapeDomain2D<T>(*this);
	}
private:
	plint cx, cy;
	plint radiusSqr;
};

然后，这个类的一个实例作为参数提供给defineDynamics:

defineDynamics( lattice, lattice.getBoundingBox(),
				new CylinderShapeDomain2D<T>(cx, cy, radius),
				new BounceBack<T,DESCRIPTOR> );

第二个参数是Box2D类型的，它可以用来提高这个函数调用的效率:在这个框和域函数指定的域之间的交集上的单元格上赋予反弹节点属性。
如果您的域被指定为一组更简单的域的并集，那么您可以在每个更简单的域上迭代地使用defineDynamics。如果域是简单域的交集或其他几何运算，则需要在域函数的定义中使用布尔运算符。

从布尔掩码返回域
在本节中，假设定义域的几何形状是由外部源规定的。例如，当您模拟多孔介质并拥有所述多孔材料的扫描数据时，就会出现这种情况。这种最简单的方法是将数据存储为ASCII文件，该文件通过0和1区分流动节点和实体节点，用空格、制表符或换行符分隔。该文件表示常规数组的内容，只存储原始数据，不存储关于矩阵维的信息。数据布局必须与Palabos相同:z坐标的增量表示文件中的连续进度(在3D中)，接着是y坐标，然后是x坐标。这些数据只是从文件中刷新到一个标量字段中。下面的代码摘自例子 examples/codesByTopic/io/loadGeometry.cpp:

MultiScalarField2D<T> boolMask(nx, ny);
plb_ifstream ifile("geometry.dat");
ifile >> boolMask;

请注意，目前没有为此类I/O操作实现错误检查。您有责任确保标量字段的大小与几何文件中的数据大小一致。下一步，使用函数defineDynamics的一个版本实例化回跳节点:

defineDynamics(lattice, boolMask, new BounceBack<T,DESCRIPTOR>, true);

这个函数目前仅使用相同的数据类型(在本例中为T)为多个块定义，这解释了为什么scalar-field布尔掩码的类型是T而不是bool，正如人们所期望的那样。标量场不是人们所期望的bool，而是真正的T类型。

解释说明

对于边界条件，操作方法都在这了，基本没什么遗漏，反弹边界条件是通过动态类设置的，是可以覆盖设置的，OnLatticeBoundaryConditionXD类，如果没看错，是通过使用枚举值调用特定的函数来实现的，不能重复设置，请注意，一般情况下，因为我主要是做多孔介质方面，一个一个边角点设置边界没有可能，我选用的就是最基本的反弹边界条件。
因为并没有做到，我并没有继续翻译后面一小段的曲面边界的内容。如果有需要，请自行观看实验。
还应注意，不同的模型对于边界的处理也是不同的，相应的调用的枚举值也是不一样的，需要各位在仔细阅读相关的模型案例后，自行判断使用方法。

完