openfoam中cyclic周期性边界的问题

李东岳

网格太多了我这跑不动，如果可以减少到1万个网格我可以跑跑看看。

余正东

@李东岳好的，感谢李老师的回复，可以简化成二维问题试一试~

xiaofenger

@余正东交流当然可以啊，不过关于cyclic边界条件我也不是特别懂。。。

余正东

@xiaofenger 你好，我也是这方面的新手，我就是想问问你的算例最后怎么跳出时间步长特别小的情况，我今天又重新跑了算例，发现还是有问题！qq：768620698，谢谢！

余正东

@xiaofenger 你好，我还想问一下请问一下你用的是cyclic还是cyclicAMI边界，另外你用的是VOF模型吗？谢谢~

余正东

@YHDTHU 请问一下这个可能是由于MULES算法的原因吗，导致在周期性边界处，不能收敛

yhdthu

@余正东不知道你在边界处有没有加密网格？

余正东

@yhdthu 在边界处局部加密吗

yhdthu

@余正东可以把网格介绍一下

余正东

@yhdthu 我的算例很简单

blocks
(
    hex (0 1 2 3 4 5 6 7) (160 320 160) simpleGrading (1 1 1)
);

程迪

cyclic BC是个大坑，因为OF的理论和编程经常不自洽。感觉网上从来没说清楚过。刚才我才搞明白。

首先，face addressing的owner, neighbour和ldu addressing的lower, upper并不天然等价。face addressing描述的的网格的几何拓扑关系。而ldu addressing需要描述的的是数值拓扑关系，如果没有耦合边界，二者是等价的，但是如果有耦合边界，就有了face addressing描述不了的off-diagonal term。这就有两种处理方式：
在OF+1706中，有两种lduAddressing，包括fvMeshLduAddressing和fvMeshPrimitiveLduAddressing。
- 实际上fvMeshLduAddressing是将face addressing和ldu addressing等价了，内部面的owner就是lower，neighbour就是upper。这也是最最常用的类型，可以称之为普通ldu。
- 而fvMeshPrimitiveLduAddressing（这里的Primitive的意思是not developed or derived from anything else ）中的lower可以添加owner之外的相互作用关系。但是这个类位于overset目录下，可见它主要用于重叠网格的情况，可以称之为overset ldu。
而coupled BC是实现在普通ldu中。依靠两个类：coupledFvPatch和coupledFvPatchField，其中关键是他们各自继承了lduInterface和lduInterfaceField，并实现了相关函数。所以其实interface和coupled BC是几乎同样的概念。
从概念上讲，对于普通LDU，在有interface的情况下虽然在lduMatrix中存了upper,diag,lower矩阵系数，但是还有额外耦合项，也就是说 $A = L + D + U + C \neq L + D + U$ ， $C$ 就是额外耦合项，由于数据结构设计上的限制，不能存到普通ldu的矩阵中，需要单独处理。
处理方式也比较简单， $C$ 因为也是代表cell-cell相互作用，所以主要是非对角的，因为对角的部分完全可以包含在 $D$ 中，而且源项部分也可以容纳在原来的源项里，可以分裂为上下两部分， $C = C_{L} + C_{U}$ 。在lduMatrix.Amul()的操作中，计算 $A \cdot x$ 的操作被分为 $D \cdot x + (L + U) \cdot x + C_{U} \cdot x_{n}$ ：
- initMatrixInterfaces: $C_{L} \cdot x_{o}$ ，其中 $x_{o}$ 是 $x$ 中属于耦合边界的内侧单元，计算结果要送到另一侧去；
- all cells: ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];: $D \cdot x$
- all faces: ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];: $L \cdot x$
- all faces: ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];: $U \cdot x$
- updateMatrixInterfaces: $C_{U} \cdot x_{n}$ ，其中 $x_{n}$ 只包含 $x$ 中属于耦合边界的外侧单元；
- 这种 $C \cdot x$ 叫做patch Contribution。 $C_{L}, C_{U}$ 叫interfaceUpper, interfaceLower。
问：为什么 $C_{L}, C_{U}$ 要分开处理？答：因为在存在processor边界时， $x_{n}$ 并不在本地，要通过通信获取结果，或者获取系数。

问： $C_{L}$ 和 $C_{U}$ 存在哪儿？答：

 class fvMatrix
:
    public refCount,
    public lduMatrix
{
//...
        //- Boundary scalar field containing pseudo-matrix coeffs
        //  for internal cells
        FieldField<Field, Type> internalCoeffs_; //注意这里不是引用，是真货
 
        //- Boundary scalar field containing pseudo-matrix coeffs
        //  for boundary cells
        FieldField<Field, Type> boundaryCoeffs_; //都初始化为0；
//...
}

coupled BC, interface, constraint BC等的关系：
interface就是coupled BC，constraint BC是OF文档中的分类，其实包含了empty, wedge,symmetry等几何约束类的BC和coupled BC。coupled BC中比较重要的两类是cyclic和processor。

总的来看cyclic边界条件本身的实现应该没有问题，还挺巧妙的。

回到这个问题，我觉得可能是CFL数计算的部分可能在cyclic BC处有bug。
Courant数计算方式来看，普通Courant数和其他一样，额外多了一个Interface Courant数。里面多了个这玩意儿：

     scalarField sumPhi
    (
        mixture.nearInterface()().primitiveField()
       *fvc::surfaceSum(mag(phi))().primitiveField()
    );
 
    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 
    meanAlphaCoNum =
        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
 
// mixture.nearInterface()()，定义是：
 
 Foam::tmp<Foam::volScalarField>
Foam::interfaceProperties::nearInterface() const
{
     return pos(alpha1_ - 0.01)*pos(0.99 - alpha1_);
}

可能是相边界速度太大了吧。估算一下有上万了，你的平均Courant这么小，最大Courant这么大。。。你又是这么均匀的网格。。

dzw05

@程迪大佬可以详细的解释一下cyclic的计算过程中吗？比如这个CL和CU到底对应哪一项，initMatrixInterfaces怎么计算了Cl*x0.

Tong

@dzw05
最近在研究processorFvPatch，来挖个坟补充下~~
对于并行耦合边界条件，在方程组迭代求解中的“耦合”过程在矩阵-向量乘法（Amul）中的实现如下：
1、initMatrixInterfaces：只负责初始化边界数据，将耦合边界内侧（本地）的数据传送到外侧（相邻进程）。
2、本地数据与向量相乘，参考程迪的描述。
3、updateMatrixInterfaces：使用第一步初始化得到的数据来耦合外侧变量（ $C \cdot x_{n}$ ）.
并行边界条件中，对于一个本地进程来说耦合过程是单向的，只将外侧数据耦合到本地，但是对于一个耦合边界面来说耦合过程是双向的，即面两侧的进程都需要耦合外侧数据（ $C = C_{L} + C_{U}$ 一侧是 $C_{L}$ 另一侧是 $C_{U}$ ）。
耦合边界类使用的接口是initInterfaceMatrixUpdate，updateInterfaceMatrix，这两个接口的实现在对应的子类边界条件中，如何想研究何以找找。

Micro

@Tong 非常期待您的后续关于处理器边界的解读，我也正在关注OF中的处理器边界，想请问一下，您知道对于CL和CU这两个系数，他们分别存放在哪里吗？以及他们的索引是存放在哪里吗，我感觉程迪老师介绍的好像有些出入，想请教一下您，谢谢！

Tong

@Micro
好像我在知乎上你的PCG的文章下面回复过，嘿嘿~
CL和CU对于耦合边界来说没有区分，都是interface上的系数，Amul()里为psi乘这个interface系数的方法是和interface绑定的，具体来说Amul()中的updateMatrixInterfaces会遍历所有的coupled interface 然后每个interface里都施行乘法操作，在doxygen里面顺着updateMatrixInterfaces一路查下去就找到了。

Micro

@Tong 是的哈，系数CL和CU都存处在interface中的boundarycoeffs中，不过我现在没搞清楚局部索引和全局索引是储存在哪里，尤其是他们是怎么相互转化的。。我一直对AX=b感兴趣，如果能彻底弄清楚OF中的矩阵元素和他们索引的相关接口，我就可以把外部的线性代数库集成到OF中。这是我很想做的。。

Tong

@Micro
之前好像记得国防科大有人用PETSc来求解openfoam线性方程组的，文章里写了怎么把OF的ldu Matrix与边界条件的架构和PETSc的线性方程组求解器融合在一起的，你可以搜来看看。我感觉对于of的系数矩阵A是按照lduMatrix的方式存的，如果外部求解器稀疏矩阵格式也是类似lduMatrix的，可以不对数据格式进行转化，否则需要将lduMatrix和转化为外部线性库支持的求解格式。我最近在做重叠网格的并行，希望可以多交流~~。

Micro

@Tong 好的！谢谢，您说的文章我也看过，不过涉及到具体细节他没有将，OF2006中已经实现了petsc2foam了。我正在看源码，在这个网址上：https://develop.openfoam.com/modules/external-solver/-/blob/develop/src/petsc4Foam/utils/petscLinearSolverContext.H
我最终的想法是实现一个自己的预处理器或者矩阵求解器然后继承到OF中。

Tong

@Micro
foam-extend 里面好像也有不少新的求解器，或许可以参考一下
https://github.com/Shadow-fax/foam-extend-4.1/tree/93fffffdb7453736d95f5535f9c8eed1054e48f6/src

最近发现对于全耦合求解的矩阵（foam-extend里的coupledMatrix）的确需要针对物理问题的新预条件器来加速收敛。

Micro

@Tong 好的谢谢。我好好看一下你提到的foam-extend，您提到：“对于全耦合求解的矩阵（foam-extend里的coupledMatrix）的确需要针对物理问题的新预条件器来加速收敛”，您的意思是预处理（相比分离求解器来讲）对于耦合求解器来说更加重要吗？

	class fvMatrix
	:
	public refCount,
	public lduMatrix
	{
	//...
	//- Boundary scalar field containing pseudo-matrix coeffs
	// for internal cells
	FieldField<Field, Type> internalCoeffs_; //注意这里不是引用，是真货

	//- Boundary scalar field containing pseudo-matrix coeffs
	// for boundary cells
	FieldField<Field, Type> boundaryCoeffs_; //都初始化为0；
	//...
	}

	scalarField sumPhi
	(
	mixture.nearInterface()().primitiveField()
	*fvc::surfaceSum(mag(phi))().primitiveField()
	);

	alphaCoNum = 0.5gMax(sumPhi/mesh.V().field())runTime.deltaTValue();

	meanAlphaCoNum =
	0.5(gSum(sumPhi)/gSum(mesh.V().field()))runTime.deltaTValue();

	// mixture.nearInterface()()，定义是：

	Foam::tmp<Foam::volScalarField>
	Foam::interfaceProperties::nearInterface() const
	{
	return pos(alpha1_ - 0.01)*pos(0.99 - alpha1_);
	}

CFD中文网

openfoam中cyclic周期性边界的问题