openfoam中cyclic周期性边界的问题
-
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矩阵系数,但是还有额外耦合项,也就是说
, 就是额外耦合项,由于数据结构设计上的限制,不能存到普通ldu的矩阵中,需要单独处理。 -
处理方式也比较简单,
因为也是代表cell-cell相互作用,所以主要是非对角的,因为对角的部分完全可以包含在 中,而且源项部分也可以容纳在原来的源项里,可以分裂为上下两部分, 。在lduMatrix.Amul()
的操作中,计算 的操作被分为 :initMatrixInterfaces
: ,其中 是 中属于耦合边界的内侧单元,计算结果要送到另一侧去;- all cells:
ApsiPtr[cell] = diagPtr[cell]*psiPtr[cell];
: - all faces:
ApsiPtr[uPtr[face]] += lowerPtr[face]*psiPtr[lPtr[face]];
: - all faces:
ApsiPtr[lPtr[face]] += upperPtr[face]*psiPtr[uPtr[face]];
: updateMatrixInterfaces
: ,其中 只包含 中属于耦合边界的外侧单元;- 这种
叫做patch Contribution
。 叫interfaceUpper
,interfaceLower
。
-
问:为什么
要分开处理?答:因为在存在processor边界时, 并不在本地,要通过通信获取结果,或者获取系数。 -
问:
和 存在哪儿?答:-
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
最近在研究processorFvPatch
,来挖个坟补充下~~
对于并行耦合边界条件,在方程组迭代求解中的“耦合”过程在矩阵-向量乘法(Amul
)中的实现如下:
1、initMatrixInterfaces
:只负责初始化边界数据,将耦合边界内侧(本地)的数据传送到外侧(相邻进程)。
2、本地数据与向量相乘,参考程迪的描述。
3、updateMatrixInterfaces
:使用第一步初始化得到的数据来耦合外侧变量( ).
并行边界条件中,对于一个本地进程来说耦合过程是单向的,只将外侧数据耦合到本地,但是对于一个耦合边界面来说耦合过程是双向的,即面两侧的进程都需要耦合外侧数据( 一侧是 另一侧是 )。
耦合边界类使用的接口是initInterfaceMatrixUpdate
,updateInterfaceMatrix
,这两个接口的实现在对应的子类边界条件中,如何想研究何以找找。 -
@Tong 好的!谢谢,您说的文章我也看过,不过涉及到具体细节他没有将,OF2006中已经实现了petsc2foam了。我正在看源码,在这个网址上:https://develop.openfoam.com/modules/external-solver/-/blob/develop/src/petsc4Foam/utils/petscLinearSolverContext.H
我最终的想法是实现一个自己的预处理器或者矩阵求解器然后继承到OF中。 -
@Micro
foam-extend 里面好像也有不少新的求解器,或许可以参考一下
https://github.com/Shadow-fax/foam-extend-4.1/tree/93fffffdb7453736d95f5535f9c8eed1054e48f6/src最近发现对于全耦合求解的矩阵(foam-extend里的coupledMatrix)的确需要针对物理问题的新预条件器来加速收敛。
-
@Micro
个人观点,线性方程组求解来说,感觉目前主要用的是克雷洛夫子空间方法+多重网格方法这两个大类,已经能在理论上比较好的解决求解方程组了(复杂度大概在O(N^2)~O(NlogN))。当然对于不同线性代数库对这些方法的实现效率有高有底(针对体系结构和并行的优化可能不同),但是感觉如果优化到位大概也就是到那个程度了,在同一套系统上不会有数量级的差距。
我之前做实验发现,对于有些多物理量耦合在一起的分块矩阵,通用的预条件器好像对降低条件数作用比较有限(主要也是由于物理上的特殊性造成的),针对某一类物理问题,似乎需要一些特殊的预条件器来有效的降低条件数,听说预条件这方面的论文也挺多的,不过比较偏数学,我不太懂。 -
@Tong 老师,您好。我想请教您一下,interface是针对周期边界、多求解域(如旋转机械)、并行边界这类复杂边界的吗,对于常规边界条件有影响吗,我理解的是边界条件的离散系数已经添加到lduMatrix中的对角线系数和源项中了。在测试中发现initMatrixInterfaces和updateMatrixInterfaces对于X和B都没有影响啊
28/32