T
hello李老师~ 最近看了看这本笔记，发现有的地方有点小问题~~~。
（1）多孔介质部分，P115页，公式7.58的以及后文的注解中，您将D解释为渗透率(permeability or k)，这应该是有些不准确的，根据达西公式描述的压力速度关系，渗透率的量纲为m^2 (https://www.calculator.org/properties/permeability.html)，通常左乘压力梯度(可以参考[1]中达西公式描述(Equation-(1)))，OF中的多孔介质模型实现是在动量方程添加阻力源项来实现的，与达西公式不同的是：在公式7.59中，D左乘了速度(U)，大小为渗透率的倒数，量纲为m^-2，为了避免误解，在这里解释D为达西阻力系数，F为 Forchheimer 阻力系数更合适，多孔介质渗透率(k)越大对应的达西阻力系数（D）就越小，反之亦然。个人观点：将” 沉降项 “翻译成”阻力源项“更好理解，也和Fluent的手册里对这种实现方法的解释一致。
(2)同时在 P125页中 7.9.6 setReference()， 对公式7.101的解释中，‘p无解’。实际上该方程的秩（rank）为 2（矩阵的第1、2行与第3行 线性相关），而未知数的个数为3，这种情况下，p有无穷多个线性相关解，这种情况下设置参考值能够将这无穷多个线性相关解退化到该参考值下的定解。同时，建议将rank翻译成’秩‘，符合国内线代教材对矩阵所在的线性空间的维数这一定义，用‘阶’的话容易和‘m*n阶’矩阵这种叫法中的‘阶’混淆造成误解。
[1]T. Wang and J. Yao, An improved embedded discrete fracture model and domain connectivity algorithms on 3D
unstructured grids, Journal of Computational Physics, 111142, doi: https://doi.org/10.1016/j.jcp.2022.111142. (求引用QAQ)
应该是精度的问题，polyMesh 的findCell 的原型为：findCell(const point & p,const cellDecomposition decompMode = CELL_TETS) 默认是使用polyMesh里的八叉树进行空间搜索的，vector每个点是一个双精度值，有效数字大约有11位置，而Info只输出了部分有效数字，感觉你直接这样写出来，就会损失精度
vector tmpKnotPoint(1.5,0.72,1.359);
vector point1(1.47727,0.694583,1.30181);
vector point2(1.5,0.694583,1.30181);
如果网格尺寸很小的话可能造成差异，
还有就是，我的记得polyMesh 网格索引树默认的精度是1e-3 如果你的网格尺寸小的话你可以调整一下，具体的你可以去看indexedOctree<>这个类的实现
@Micro
个人观点，线性方程组求解来说，感觉目前主要用的是克雷洛夫子空间方法+多重网格方法这两个大类，已经能在理论上比较好的解决求解方程组了(复杂度大概在O(N^2)~O(NlogN))。当然对于不同线性代数库对这些方法的实现效率有高有底（针对体系结构和并行的优化可能不同），但是感觉如果优化到位大概也就是到那个程度了，在同一套系统上不会有数量级的差距。
我之前做实验发现，对于有些多物理量耦合在一起的分块矩阵，通用的预条件器好像对降低条件数作用比较有限（主要也是由于物理上的特殊性造成的），针对某一类物理问题，似乎需要一些特殊的预条件器来有效的降低条件数，听说预条件这方面的论文也挺多的，不过比较偏数学，我不太懂。
@Micro
foam-extend 里面好像也有不少新的求解器，或许可以参考一下
https://github.com/Shadow-fax/foam-extend-4.1/tree/93fffffdb7453736d95f5535f9c8eed1054e48f6/src
最近发现对于全耦合求解的矩阵（foam-extend里的coupledMatrix）的确需要针对物理问题的新预条件器来加速收敛。
@dzw05
最近在研究processorFvPatch，来挖个坟补充下~~
对于并行耦合边界条件，在方程组迭代求解中的“耦合”过程在矩阵-向量乘法（Amul）中的实现如下：
1、initMatrixInterfaces：只负责初始化边界数据，将耦合边界内侧（本地）的数据传送到外侧（相邻进程）。
2、本地数据与向量相乘，参考程迪的描述。
3、updateMatrixInterfaces：使用第一步初始化得到的数据来耦合外侧变量（$C\cdot x_n $）.
并行边界条件中，对于一个本地进程来说耦合过程是单向的，只将外侧数据耦合到本地，但是对于一个耦合边界面来说耦合过程是双向的，即面两侧的进程都需要耦合外侧数据（$C=C_L+C_U$ 一侧是$C_L$另一侧是$C_U$）。
耦合边界类使用的接口是initInterfaceMatrixUpdate，updateInterfaceMatrix，这两个接口的实现在对应的子类边界条件中，如何想研究何以找找。
今天发现这个问题折腾了小半天也没解决，恳请大家帮忙分析一下：
非稳态计算过程中需要连续的向模型中注入流体，由于尺度的问题无法使用“显式”的入口，而是使用一个“隐式”的方法：添加一个体积源项到压力方程中实现注入；
体积源项定义：
，其实就是想在某些网格的连续方程上添加一个源项，来实现流体注入:
volScalarField wellSource(
IOobject("wellSource",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE),
mesh,
dimensionedScalar("Scalar", dimVol / dimTime, Zero)
);
需要添加体积源项的网格对应的wellSource有数值，否则为零，注入量是一个固定不变的值。
整体求解流程类似PISOFoam 前后代码就没有列出来，压力方程源项添加如下：
fvScalarMatrix pEqn((fvm::laplacian(rTU, p) == fvc::div(phiHbyA))-(wellSource/mesh.V()));
但是在计算过程出现问题是：每一个时间步求解一次的压力方程后压力都会增加，对应的库朗数也一直增加（如下）：
Courant Number mean: 0.000925584 max: 0.19292
Courant Number mean: 0.00103634 max: 0.331209
Courant Number mean: 0.00132796 max: 3.74959
Courant Number mean: 0.00177495 max: 4.40821
Courant Number mean: 0.00422436 max: 24.4105
Courant Number mean: 0.0173964 max: 47.9946
Courant Number mean: 0.0974709 max: 287.174
Courant Number mean: 0.679524 max: 1792.34
Courant Number mean: 39.8774 max: 147650
...........
Courant Number mean: 8.97372e+11 max: 6.71492e+15
看上去是每一个时间步的源项和上一个时间步的结果都在累加，从压力方程上每个时间步计算的过程来看的确也是这样的。
从物理过程的角度来说，连续的向模型中注入流体，每个时间步的压力方程都应该有源项存在的，如果源项大小不变，压力应该是稳态的，
不知道要如何避免这种压力不断增加的结果~,希望大家提提意见。