首先，对于压力方程，因为其只有 laplacian项，所以可以推导出来这样的压力方程的系数满足 “网格单元对应的对角矩阵系数等于其余影响的单元系数之负和” 这个特点。
其次，我们来看 setReference 这个函数做了什么。查看 fvMatrix 类中 setReference 函数的实现，如下：
template<class Type>
void Foam::fvMatrix<Type>::setReference
(
const label celli,
const Type& value,
const bool forceReference
)
{
if ((forceReference || psi_.needReference()) && celli >= 0)
{
source()[celli] += diag()[celli]*value;
diag()[celli] += diag()[celli];
}
}
可以看到，这个函数确实是把对应网格的对角元素值翻倍了。 setReference 函数的目的是把指定网格（对应参数 celli）的值设置为指定的值（对应参数 value），可以为什么这样实现就可以达到这个目的呢？原理如下：
在不可压缩流体计算中，我们求解的是压力泊松方程（Pressure Poisson Equation），它通常形如：
$$ \nabla \cdot (\frac{1}{A_p} \nabla p) = \nabla \cdot (\frac{H(\mathbf{U})}{A_p}) $$
这里 p 是压力，A_p 是速度方程离散后中心网格的系数，H(U) 是一个与速度相关的项。
当我们将这个方程离散化后，会得到一个大型的线性方程组：
$$ \mathbf{A} \mathbf{p} = \mathbf{b} $$
其中 A 是系数矩阵，p 是我们要求解的、包含所有网格压力的向量，b 是源项向量。
关键问题：如果计算域的所有边界都是第二类边界条件（Neumann BC），比如 zeroGradient（法向梯度为零），这意味着我们只规定了压力的梯度，而没有规定压力的绝对值。
从物理上讲，不可压缩流的压力是相对的。如果 p 是一个有效的解，那么 p + C（其中 C 是任意常数）也同样是一个有效的解，因为它不会改变压力的梯度（∇(p+C) = ∇p）。
因此，我们必须通过指定一个点的压力值来消除这种不确定性，为整个压力场提供一个“锚点”或“基准”。
setReference 的目标是强制让某个特定网格 celli 的压力 p[celli] 在求解后等于我们设定的参考值 value。
也就是说，我们想让方程组的解满足：
$$ p_{\text{celli}} = \text{value} $$
3. 代码实现与数学原理剖析
我们再来看这两行代码，并结合线性方程组 Ap=b 在第 celli 行的具体形式来分析。
对于网格 celli，其对应的线性方程（未修改前）是：
$$ A_{ii} p_i + \sum_{j \in N(i)} A_{ij} p_j = b_i $$
其中：
i就是cellip_i是网格i的压力A_{ii}是对角线元素，对应 OpenFOAM 代码中的diag()[celli]A_{ij}是非对角线元素，表示网格i与其邻居网格j的相互影响b_i是源项，对应source()[celli]
现在我们来分析这两行代码对这个方程做了什么修改：
核心代码:
source()[celli] += diag()[celli]*value;
diag()[celli] += diag()[celli];
第一步：source()[celli] += diag()[celli]*value;
这行代码修改了源项 b_i。修改后的新源项 b'_i 为：
$$ b_i^{'} = b_i + A_{ii} \cdot \text{value} $$
此时，方程变为：
$$ A_{ii} p_i + \sum_{j \in N(i)} A_{ij} p_j = b_i + A_{ii} \cdot \text{value} $$
第二步：diag()[celli] += diag()[celli];
这行代码修改了对角线元素 A_{ii}。修改后的新对角线元素 A'_{ii} 为：
$$ A_{ii}^{\prime} = A_{ii} + A_{ii} = 2A_{ii} $$
现在，第 i 行的方程最终变成了：
$$ (2A_{ii}) p_i + \sum_{j \in N(i)} A_{ij} p_j = b_i + A_{ii} \cdot \text{value} $$
让我们来整理一下这个修改后的方程。我们可以把 (2A_{ii}) p_i 拆开：
$$ (A_{ii} p_i + \sum_{j \in N(i)} A_{ij} p_j) + A_{ii} p_i = b_i + A_{ii} \cdot \text{value} $$
也就是说，当调用了 setReference 函数后，求解的线性方程组发生了变化。通过求解这个修改后的方程，能让方程趋于收敛到 $p_i = \text{value}$ 这个解，此时
$$ A_{ii} p_i = A_{ii} \cdot \text{value} $$
因此也会使得修改前的
$$ A_{ii} p_i + \sum_{j \in N(i)} A_{ij} p_j \approx b_i $$
也同时成立。
通过这样就实现了让指定的网格的解等于指定值，同时又不违反从 Laplacian 项离散后得到的方程。
